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ABSTRACT 


Algorithms  recommended  for  inclusion  in  an  update  of 
the  AD  real-time  range  safety  control  program,  P2457,  are 
developed.  The  developed  algorithms  include  ones  for  large 
and  small  beta  fragment  impact  prediction  by  table  look-up 
procedures,  adaptive  exponential  digital  filtering  of  radar 
data,  and  the  optimal  combining  of  radar  and  telemetered 
navigation  system/altimeter  system  data  for  use  in  impact 
prediction.  The  algorithms  are  presented  in  a  form  conducive 
to  fast,  efficient  incorporation  in  the  P2457  update. 
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INTRODUCTION 


During  Task  SER  82-1,  algorithms  which  should  be 
included  in  an  updated  version  of  the  real-time  range  safety 
control  program  (P2457)  for  the  VAX  11/780  real-time  computer 
configuration  were  identified.  The  identified  algorithms  were 
of  three  types: 

a)  Algorithms  being  utilized  in  P2457,  requiring 
no  modifications, 

b)  Algorithms  utilized  in  P2457,  requiring  modifi¬ 
cations,  and 

c)  New  algorithms  not  currently  included  in  P2457. 

Among  the  new  algorithms  identified  and  recommended  were  the 
following: 

1)  Because  vacuum  instantaneous  impact  predictions 
( I  IPs )  were  found  to  be  inappropriate  for  most 
tactical  weapon  tests  and  the  current  P2457 
drag/wind  IIP  algorithm  was  found  to  be  too  slow 
for  multiple  weapon  tests,  the  use  of  table 
lookup  IIP  algorithms  for  the  largest  and  smallest 
beta  fragments  of  a  destroyed  vehicle  were 
recommended  for  inclusion  in  a  P2457  update. 

2)  Because  the  current  radar  digital  filter  algorithm 
used  in  P2457  was  found  to  have  several  deficiencies, 
including  instability,  a  fast,  stable  exponential 
weighted  filter  algorithm  was  recommended. 
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3)  Because  of  radar  data  degradation  at  the  low 
elevation  angles  common  on  advanced  AD  weapon 
tests,  the  inclusion  of  algorithms  for  processing 
telemetered  inertial  navigation  and  altimeter 
data  was  recommended. 

The  purpose  of  this  study  was  to  develop  the  above 
three  algorithms  and  present  them  in  a  form  conducive  to  fast, 
efficient  incorporation  into  a  P  2  4  5  7  update  on  the  VAX  11/780 
configuration.  Other  required  algorithms  identified  in  Task 
SER  82-1  will  be  developed  later. 
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2.0 


SUMMARY  AND  RECOMMENDATIONS 


2.1  Summary 

The  following  algorithms  were  developed  during  this 
study  and  presented  in  a  form  suitable  for  efficient  incorporation 
in  a  P2457  update  on  the  AD  VAX  11/780  real-time  computer 
conf i gura ti on  . 

1)  Table  lookup  algorithms  for  generating  the  large 

and  small  beta  fragment  vertices  of  a  debris  triangle. 

2)  A  stable,  fast,  adaptive  exponential  filter  algorithm 
to  replace  the  current  unstable  radar  data  filter 
algorithm  of  P2457. 

3)  Algorithms  for  combining  radar  data  with  telemetered 
INS  and  altimeter  data  to  obtain  smooth,  unbiased 
trajectory  state  vector  estimates. 

The  above  developed  algorithms  can  be  included  in  the 
development  of  a  type  B5  computer  program  specification,  as 
recommended  in  Task  SER  82-1. 

2 . 2  Recommendati ons 

It  is  recommended  that  the  impact  prediction,  adaptive 
exponential  filtering,  and  radar/TM  data  combinative  algorithms 
developed  during  this  study  be  incorporated  in  the  planned 
update  of  the  real-time  range  safety  control  program,  P2457.  The 
algorithms  should  be  included  in  a  type  B5  computer  program 
specification  for  the  P2457  update. 
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3.0 


REAL-TIME  IMPACT  PREDICTION  ALGORITHMS 


3 . 1  Background 

Several  reasons  why  the  current  impact  prediction 
algorithms  in  P2457  should  be  replaced  were  cited  in  Reference 
1.  The  primary  reason  given  for  discarding  the  vacuum  impact 
predictor  was  that  its  predictions  have  become  increasingly 
unrealistic  and  too  restrictive  for  advanced  tactical  weapon 
tests,  which  are  conducted  primarily  in  the  lower  regions  of 
the  atmosphere.  For  example,  consider  the  following  case, 
which  is  highly  possible  for  near-future  weapons: 


•r-  -o 
C  C 

►— *  o 
<_> 


A 1 1 T tude 
Speed 

Flight  path  angle 

Fragment  ballistic  coe f f i c i en t ( g ) 

True  down  range  impact  distance 
Vacuum pred  i  cted  down  range  distance 


5000  feet 
Mach  2.7 

30°  above  horizontal 
500 

36,500  feet 
250 , 000  feet 


It  can  easily  be  seen  that  use  of  the  vacuum  impact  prediction 
to  determine  test  boundaries  would  be  unnecessarily  restrictive. 

The  primary  reason  given  in  Reference  1  for  discarding 
the  current  drag  and  wind  corrected  impact  predictor  algorithm 
was  that  the  required  numerical  integration  is  too  slow  to  be 
used  for  predicting  impacts  on  tests  involving  multiple  weapons 
and  multiple  fragments.  Each  prediction  can  consume  as  much  as 
200  milliseconds  on  the  CDC  6600,  and  the  time  required  could  be 
greatly  increased  on  the  VAX  11/780,  where  the  use  of  double 
precision  might  also  be  necessary. 
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Reference  1  recommended  that,  for  an  update  of  P2457 
on  the  VAX  11/780  configuration,  table  look-up  algorithms 
similar  to  those  proposed  in  Reference  2  should  be  developed. 

Two  types  of  impact  prediction  algorithms  are  needed:  a)  an 
algorithm  for  predicting  the  impact  point  of  the  fragment 
having  the  largest  6  (or  the  intact  vehicle)  and  b)  an  algorithm 
for  predicting  the  impact  point  of  the  smallest  6  fragment  which 
is  hazardous.  It  was  also  recommended  in  Reference  1  that  the 
two  predicted  impact  points  be  presented  on  a  graphics  display 
as  two  vertices  of  a  debris  triangle,  as  shown  in  Figure  3.1.1. 

One  of  the  purposes  of  this  study  is  to  develop  the  large  and 
small  6  impact  prediction  algorithms  in  a  form  conducive  to 
inclusion  in  an  update  of  P2457  on  the  VAC  11/780  configuration. 

3.2  Lagrangian  Interpolation  and  Extrapolation  of 

Multivariate  Tables 

Table  look-up  impact  predictors  require  the  use  of 
interpolation  algorithms.  In  developing  an  interpolation 
algorithm,  one  must  consider  the  following: 

1)  Number  of  independent  variables 

2)  Spacing  of  independent  variables 

3)  Degree  of  interpolation  polynominal 

4)  Need  for  extrapolation  outside  the  table 

Based  upon  Reference  2,  it  is  concluded  here  that  the  required 
interpolation  algorithm  must  be  able  to  accept  three  unevenly  spaced 
independent  variables  (trivariate  interpolation).  Selection  of 
the  interpolation  polynominal  degree  involves  a  tradeoff  between 
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using  linear  interpolation  of  a  table  requiring  a  large 
block  of  storage  and  higher  degree  interpolation  requiring 
less  storage.  In  order  to  make  optimal  use  of  the  8K  byte  VAX 
11/780  cache  memory,  the  use  of  small  tables  would  be  desirable. 
Lagrangian  interpolation  is  a  widely  used  method  for  performing 
polynominal  interpolation  and  extrapolation  of  tables  having 
unevenly  spaced  variables.  For  second  degree  polynominal 
interpolation  and  extrapol ati on ,  the  Lagrange  formula  is: 

f(X)  =((X-Xi)(X-X2)  /  (Xo-Xi ) (X0-X2) )  f(Xo) 

+  ((X-X0)(X-X2)  /  (X1-X0) (X1-X2 ) )  f ( X 1 )  (3.2.1) 

+  (( X - X 0 ) (X-Xi )  /  (X2-X0) (X2-X1 ) )  f(X2) 

where 

X  =  coordinate  of  desired  interpolated  point 
f(X)  =  the  desired  interpolated  function  approximation 
Xc,Xi,X2  =  three  values  of  the  independent  variable 
f ( X 0 ) , f ( X  1 )  ,  f ( X2  )  =  functional  values  for  X  =  Xo,Xi,X2 

Available  bivariate  and  univariate  second  degree 
Lagrangian  interpolation  subroutines  (References  12  and  13)  were 
combined  to  form  the  required  trivariate  subroutine.  A  listing 
of  the  subroutine  (LAGRAN)  is  provided  in  Appendix  A.  Since  a 
primary  reason  for  recommending  table  look-up  drag/wind  corrected 
impact  prediction  was  computational  speed,  a  large  three  dimen¬ 
sional  table  (6x8x10=480  elements)  was  input  to  LAGRAN  on  the 
CDC  6600  computer  and  the  interpolation  was  timed.  It  was  found 
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that  100  interpolations  could  be  performed  in  0.045  seconds, 
or  in  less  than  0.5  milliseconds  per  interpolation.  Thus,  the 
table  look-up  algorithm  is  as  much  as  400  times  as  fast  as  the 
current  real-time  drag/wind  corrected  impact  predictor  (Sub¬ 
routine  IMPACT)  and  should  be  suitable  for  use  on  tests  re¬ 
quiring  predictions  for  several  objects. 

3.3  T rajectory  Integrator 

The  current  drag/wind  corrected  impact  prediction 
subroutine  (IMPACT)  was  evaluated  in  Reference  12  and  was  found 
to  provide  satisfactory  impact  predictions  for  tactical  weapon 
tests,  although  minor  modifications  were  recommended.  A  modified 
version  of  IMPACT  and  its  auxiliary  subroutines  are  recommended 
for  inclusion  in  the  P2457  update;  however,  they  will  be  used 
to  develop  the  required  large  and  small  B  impact  prediction 
tables  rather  than  to  predict  impacts  during  real-time  execution. 
For  the  development  of  the  tables  the  following  modifications  are 
recommended : 

1)  The  hazardous  small  6  body  requiring  impact 

predictions  on  some  tests  (e.g.,  MGD,  FIREB0LT) 
will  be  the  parachute  and  its  load.  Because 
of  the  slow  rate  of  descent,  the  integration 
of  parachute  trajectories  to  impact,  using  the 
IMPACT  subroutine,  is  extremely  time  consuming. 

A  version  of  the  fast  parachute  prediction 
algorithm  presented  in  Appendix  B  should  be  added 
to  the  IMPACT  codes  and  be  used  for  generating 
parachute  impact  prediction  tables. 
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2 )  The  current  version  of  IMPACT  utilizes  an 

exponential  atmospheric  density  model  (recently 
modified  from  a  quadratic  model).  The  densities 
given  by  the  exponential  model  can  differ  from 
the  true  densities  by  several  percent  and  can  also 
result  in  impact  prediction  errors  of  several 
percent.  It  is  recommended  that  the  subroutine 
be  modified  to  interpolate  atmospheric  densities 
from  the  latest  weather  deck.  The  linear  inter¬ 
polation  subroutine  INPOL,  listed  in  Appendix  B, 
can  be  used  to  perform  the  required  interpolation. 

3 . 4  Vehicle  and  Fragment  Weight,  Area  and  Drag  Data 

Vehicle  or  fragment  area  and  weight  and  coefficients 
of  drag  versus  Mach  number  are  currently  input  to  P2457  by  card 
prior  to  a  test.  For  the  P2457  update,  it  is  recommended  that 
a  library  of  this  information  be  stored  on  a  tape  or  disk  file 
and  be  selected  automatically  as  a  function  of  test  vehicle  type. 
One  of  the  recommendations  of  Task  SER  82-6  is  that  a  menu  defining 
the  test  vehicles  be  presented  on  a  CRT  prior  to  a  test  and  be 
selectable  by  an  interactive  device  such  as  a  joy stick  controlling 
a  cursor.  Among  the  characteristics  defining  the  vehicle  are  the 
following: 

•  Vehicle  name  (e.g.,  Bomarc) 

•  Shooter/target/missile  relationships 

•  Radar  assignments 

•  Destruct  key 

If  the  82-6  recommendations  are  accepted,  it  will  be 
convenient  to  identify  the  proper  disk  file  for  the  vehicle  through 
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use  of  the  selected  name.  For  the  purpose  of  generating  the 
impact  prediction  tables,  it  will  also  be  convenient  to  in¬ 
clude  on  the  data  files  default  values  defining  the  vehicle 
pei  Tui mdnce  envelopes.  In  order  to  handle  sudden  new  require¬ 
ments,  the  P2457  update  should  also  offer  the  option  of  overriding 
the  values  on  the  library  file  through  the  keyboard. 

3.5  S m a  1 1  g  Impact  Prediction  T able s 

The  smallest  hazardous  6  fragment  impact  prediction 
tables  recommended  for  a  P2457  update  are  based  upon  the  recom¬ 
mendations  of  Reference  4.  Real-time  implementation  of  the 
algorithm  requires  the  following  steps: 

1)  Prior  to  the  test,  enter  the  proposed  modified 
version  of  IMPACT  with  the  following  data: 


a) 

b) 

c) 

d) 

e) 

f) 

g) 

h) 


Test  origin  latitude,  longitude,  height 
Flight  azimuth  =  0  0 

Cq  for  small  $  fragment  (or  reefed  and  open 
parachute  tables) 

Fragment  weight(W) 

Fragment  area  (A) 

Wind  velocity  components,  speed  of  sound, 

and  atmospheric  density  versus  altitude 

X  =Z  =VX  = V Y  =  VZ  =0 
0  0  0  0  0 

Series  of  altitudes,  hQi  i  =  l,n 
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2)  From  output  of  IMPACT,  form  the  impact  prediction 


tables: 

Altitude 

X  ( North 

Z ( Ea  s  t 

( f  eet ) 

(feet) 

( feet ) 

hl 

aX2 

AZ1 

h2 

AX2 

L> 

r^4 

ro 

h  _  AX  AZ 

n  n  n 

Here,  hn  should  be  higher  than  the  maximum  expected  altitude  during 
the  test. 


3)  During  real-time  execution,  linearly  interpolate  the 
AX,  AZ  tables  to  the  current  altitude,  h.  Subroutine 
INPOL,  Appendix  B,  can  be  used  for  this  purpose. 

4)  Approximate  the  correction  for  the  small  8  fragment 
initial  crosswind  velocity: 

AXc=  ARc  COS  (*) 

AZc=  ARc  SIN  (4-) 

*=  ARCTAN  (ZQ/X0) 

aRc=  ( VT(h) )2C0S(y)  (l+(SIN(y))/2+Ln(V0/VT(h)))/g 


VT(h)  = 


g= 


terminal  velocity  at  altitude  h  (see  Appendix 

2+i„2*z„2)1/2 

0  0  0 

ARCSIN(Yq/Vo) 

current  velocity  vector  components 

2 

acceleration  of  gravity  (32.2  ft/ sec  ) 


B) 
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5)  Sum  the  corrections  to  obtain  the  total  corrections: 


' x  I  p 

AX 

AZIp  = 

•z 

6)  Approximate  AYJp:  AYjp  =  -( ( AX  j  p ) ^  +  ( AZ j p ) 2 ) / ( 4 . 18x1 0^ ) 

7)  Enter  subroutine  XYZLLH  with  the  present  position 
latitude  and  longitude  and  with  AXjp,  AYJp,  A^jp  t0 

obtain  the  small  B  impact  prediction. 

A  macro  flow  chart  of  the  the  proposed  small  B  impact 
prediction  process  is  presented  in  Figure  3.5.1. 
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Macro  Flow  Chart  of  Small  /3  Impact  Prediction  Process 

Figure  3.5.1 
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3.6  Large  (i  fragment  Impact  Prediction 

Ideally,  the  large  (1  downranqe,  crossrange  impact 
lookup  tables  should  be  based  upon  four  independent  fragment 
initial  condition  variables: 
h  -  altitude 

u  -  flight  path  angle,  ,  ,  i  :  ' v r  above  horizontal 

V  -  velocity  magnitude,  or  speed 
vp  -  clockwise  heading  from  north 
Using  the  current  estimates  of  the  four  parameters  from  radar 
data  as  inputs,  the  four  dimensional  tables  could  be  entered 
in  real-time  to  obtain  interpolated  values  of  the  downrange  and 
crossrange  impact  distances. 

The  inclusion  of  T  as  one  of  the  independent  variables 
is  due  mainly  to  the  questionable  need  to  account  for  wind 
forces,  although  small  coriolis  force  effects  also  vary  with 
the  initial  heading.  In  the  Reference  3  study  of  large  (i  impact 
predictors,  the  effects  of  wind  on  the  predictions  was  considered 
negligible  relative  to  the  large  differences  between  the  drag- 
corrected  predictions  and  the  vacuum  predictions,  and  ^  was 
therefore  not  included  as  one  of  the  independent  table  variables. 
Analyses  conducted  during  this  study  showed  that  the  inclusion 
of  Jj  would  require,  for  some  tests,  tables  exceeding  1000 
storage  locations.  This  was  considered  excessive  relative  to 
the  impact  prediction  improvement  which  could  be  achieved; 
therefore,  it  is  one  of  the  recommendations  of  this  study  that 
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the  large  /3  tables  should  be  functions  of  only  altitude,  flight 
path  angle  and  velocity  magnitude.  For  tests  in  which  the 
nominal  trajectory  heading  is  known,  most  of  the  error  due  to 
wind  can  still  be  accounted  for  by  entering  this  nominal  heading 
into  the  impact  prediction  subroutine  before  generating  the 
tables.  Table  3.6.1  shows  some  "worst  case"  effects  of  ignoring 
the  wind  in  performing  large  (3  impact  predictions.  The  cross¬ 
range  components  are  due  mainly  to  wind.  The  tables  also  show 
that  the  vacuum  versus  drag  corrected  impact  differences  are 
much  larger  than  the  errors  due  to  ignoring  wind.  In  generating 
the  tables  with  the  IMPACT  subroutine,  a  100  ft/sec  wind  speed 
from  270°  (due  west)  was  imposed. 

The  steps  required  in  the  recommended  large  (3  impact 
prediction  procedure  are  as  follows: 

1.  A  record  containing  the  table  independent  initial 

condition  variables  for  the  test  vehicle  of  interest 
must  be  on  an  input  file.  Also  on  the  record 
should  be  the  large  /3  fragment  (or  intact  vehicle) 
weight,  effective  drag  area  and  drag  versus  mach 
number  table.  The  record  should  be  identifed 
by  test  vehicle  type  (e.g.,  Bomarc).  Based  upon 
tests  performed  during  this  study,  the  following 
numbers  of  independent  variables  of  each  category 
should  be  included  in  the  record: 
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Large/3  Fragment  Impact  Predictions 
(Wind  Profile:  Constant  100  ft/sec  from  West) 


Altitude  (feet)  -  7  values 

Flight  path  angle  (degrees)  -  7  values 

Velocity  magnitude  (ft/sec)  -  5  values 

This  will  result  in  a  generated  table  containing 

245  dependent  variables  plus  the  19  independent 

variables.  Example  of  independent  variables 

which  might  be  required  for  Bomarc: 


h 

a 

V 

500 

-10 

250 

5000 

0 

750 

15000 

10 

1000 

30000 

30 

2000 

50000 

50 

3000 

65000 

70 

75000 

90 

2.  If  wind  corrections  are  desired,  the  nominal  mean 
heading,  ij;  ,  must  be  entered  from  the  keyboard 
for  each  vehicle  requiring  impact  prediction. 

3.  Prior  to  the  test,  a  command  to  generate  the 
large  impact  prediction  tables  is  made.  The 
proper  vehicle  data  will  be  read  from  the  vehicle 
data  file  and  all  combinations  of  the  independent 
variables  will  be  used  as  inputs  to  the  modified 
IMPACT  subroutine.  If  \p  has  not  been  entered 
(Step  2  above),  wind  speed  will  be  set  to  0.0  and 
only  one  table  containing  the  downrange  distances 
as  the  dependent  variables  will  be  generated 
(default  <p  =  0°).  If  ip  is  entered,  the  available 
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wind  data  will  be  used  and  an  additional  table 
containing  the  crossrange  values  will  be  generated. 
The  tables  must  be  stored  in  arrays  compatible 
with  the  LAGRAN  subroutine  calling  sequence  (see 
Appendix  A ) . 

4.  During  real-time  execution,  after  the  vehicle 
present  latitude  (<)>),  longitude  ( X )  ,  altitude  (h), 
flight  path  angle  (c»),  heading  ( ^ )  and  velocity 
magnitude  (V)  have  been  determined  from  radar 

(or  other)  data,  subroutine  LAGRAN  will  be  called 
for  interpolation  of  the  downrange  and  crossrange 
impact  prediction  tables  (see  Appendix  A).  The 
outputs  will  be  X  (downrange)  and  Z  (crossrange). 

5.  Subroutine  XYZLLH  must  be  entered  with  the  following 
data  to  obtain  the  desired  large  /?  impact  point: 

Present  latitude  ( ) 

Present  longitude  ( X ) 

Present  altitude  (h) 

Present  heading  (<J/) 

Interpolated  X  j  p 

Interpolated  Zjp 

Computed  Yjp  =-(h  +  (Xjp  +  z|p )/ (4 . 18  x  107) 

KODE  =  3 

XYZLLH  will  return  4>ip,  xjp  and  hjp 

A  macro  flow  chart  of  the  large /3  fragment  impact 
prediction  process  is  presented  in  Figure  3.6.1. 
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Macro  Flow  Diagram  of  the  Proposed 
Large  fi  Impact  Prediction  Process 

Figure  3.6.1 
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4.0 


REAL-TIME  RADAR  DIGITAL  FILTER  ALGORITHM 


4.1  6 a ck  g  ro u n_d 

It  was  recommended  in  Reference  1  that  the  current 
real-time  radar  data  digital  filter  algorithm,  which  is  con¬ 
tained  in  subroutine  TARGET  of  P2457,  should  be  replaced  with 
a  filter  of  the  e xpo ne n t i a  1 1 y  weighted  type  if  P2457  is  updated 
on  the  VAX  11/780.  Reasons  given  for  this  recommendation  were 
the  following: 

1.  The  TARGET  filter  is  unstable,  being  completely 
unable  to  eliminate  past  errors  in  its  output 
(e.g.,  errors  due  to  roundoff).  Divergence  of 
the  filter  due  to  roundoff  errors  when  run  on  a 
32  bit  computer  is  illustrated  in  Figure  4.1.1. 

The  stable  output  of  an  exponential  filter,  using 
the  same  input  data,  is  shown  in  Figure  4.1.2. 

2.  The  TARGET  filter  requires  a  lengthy  time  period 
( norma lly  five  seconds)  to  become  reinitialized 
after  data  dropouts. 

3.  In  order  to  avoid  divergence,  the  TARGET  algorithm 
would  have  to  be  written  in  double  precision  for 
the  VAX  11/780.  The  exponential  algorithms  would 
require  only  single  precision. 

4.  As  shown  in  Reference  9,  the  TARGET  algorithm 
requires  several  times  as  many  computational  steps 
and  several  times  the  storage  of  an  exponentially 
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weighted  filter.  This  disadvantage  would  be 
increased  (by  a  factor  as  large  as  four)  because 
of  the  need  for  double  precision  on  the  VAX  11/780. 

5.  The  proposed  expansion  of  P2457  to  handle  up  to 
eight  radar  inputs  will  further  increase  the  dis¬ 
advantages  of  using  the  TARGET  algorithm,  in  terms 
of  speed  and  storage. 

The  recommendati on  of  an  exponential  filter  algorithm 
in  place  of  the  TARGET  algorithm  was  based  largely  upon  the 
results  of  two  previous  studies  (References  8  and  9).  In 
Reference  8,  the  characteristics  of  real-time  filters  which 
have  been  used  in  range  safety  work  at  the  various  test  ranges 
was  discussed.  All  of  these  filters,  in  various  ways,  attempt 
to  compromise  between 

(a)  reducing  the  ratio  of  output  noise  to 
input  noise  to  a  minimum  and 

(b)  reducing  the  bias  in  the  filter  output  to 
a  minimum. 

The  method  of  least  squares  is  a  universal  approach  to 
achieving  (a).  Editing  schemes  involving  trend  detection  with 
subsequent  modification  of  filter  weights  are  used  to  achieve 
(b).  Two  other  important  choices  often  enter  into  the  design 
of  a  range  safety  digital  radar  filter: 

1.  A  choice  of  filtering  in  the  radar's  natural 
coordinates  (azimuth,  elevation,  range  in  the 
case  of  the  FPS-16)  or  in  a  rectangular  cartesian 
coordinate  system. 
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2 .  A  choice  of  constraining  the  output  of  the 
filter  to  satisfy  a  physically  realistic  set 
of  differential  equations  (including  the  forces 
of  drag,  lift,  thrust,  gravitation,  etc.)  or 
constraining  the  output  to  follow  a  low  degree 
polynominal  . 

The  choice  of  the  real-time  filters  in  current  range 
safety  use  at  the  various  ranges  has  been  dictated  largely  by 
computer  storage  and  speed  requi rements .  As  a  result,  recursive 
exponential  least  squares  filtering  of  cartesian  coordinate  data 
with  polynominal  constraints  has  been  the  leading  choice.  This 
type  of  filtering  also  goes  by  the  names  "  a,  (3,7  filtering", 
"discounted  least  squares",  and  "fading  memory  polynominal 
filtering".  It  has  the  following  features. 

1.  It  is  fast.  For  example,  the  recursive  algorithm 
for  a  quadratic  exponential  filter  can  be  expressed 

35  \ 
c  =  X°  -  X  -  XTj  -  XT 2 

X  =  X  +  « e  ^  (4.1.1) 

X=X+XTj+8c 

X'  =  X  +  y  c 

where 

X°  is  the  current  observation, 

T^  is  the  time  increment  between  observations. 


a,  8,  y  are  the  filter  weights,  and 

X,  X,  X  is  the  current  state  vector  estimate. 
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2  .  It  requires  a  minimum  of  computer  storage. 

Only  10  memory  locations  are  required  for  the 
above  algorithm. 

3.  It  is  readily  adaptable.  The  weights,  a  ,  /3, 
and  y  can  be  expressed  as  simple  functions  of 
one  constant  O-  Using  trend  checks  on  the 
residuals  (the  e's  in  the  above  recursion)  as  a 
basis,  the  filter's  effective  memory  may  be 
easily  contracted  or  expanded  to  arrive  at  the 
desired  compromise  between  noise  reduction  and  bias 
and  stability  control. 

4.  It  gives  the  most  weight  to  the  most  recent  data. 
This  allows  quick  response  to  changes  in  target 
motion. 

5.  It  is  stable.  Due  to  the  fact  that  weight  given 
past  data  and  past  state  vector  estimates  are  a 
function  of  age,  past  input  and  computational  errors 
gradually  disappear. 

Reference  9  contained  a  comparison  of  various  test 
range  filters  in  terms  of  speed  and  storage  requirements.  The 
speed  comparisons  were  in  terms  of  the  number  of  arithmetical 
operations  required  to  fit  one  channel  of  data  (e.g.,  the  X 
position  coordinate)  to  a  second  degree  polynominal.  The  com¬ 
parisons  for  state  vector  computation  only  (no  editing,  no 
variance  estimation),  were  as  follows: 
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Algorithm 


Number  of 
Operations 


Moving  Arc  (51  points  306+ 

Moving  Arc  (N  points)  6N+ 

Current  AD  TARGET  (after  initialization)  33 

Modified  AD  TARGET  (from  Ref. 8  -  after  27 

initialization) 

Kalman  257 

Exponential  13 


The  numerical  storage  requirement  comparisons,  again 
assuming  estimation  of  position,  velocity  and  acceleration  with 
a  second  degree  polynominal  fit  to  one  input  channel,  were  as 
follows: 


A1 gori tbm 

Number  of 
Storage  Cells 

Moving  Arc  (51  points 

207 

Moving  Arc  (N  points) 

4N  +  3 

Current  AD  TARGET  (51  points  -  after 

initialization) 

127 

Modified  AD  TARGET  (51  points  -  after 

initialization) 

63 

Kalman 

64 

Exponential 

10 

Reference  9  mentioned  two  deficiencies  of  exponential 

filters: 

•  The  filters  are  not  self-starting.  They  need 
approximate  initial  state  vectors  computed  by 
some  other  method. 
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t  Overshoot  and  undershoot  settling  problems 

are  common  after  sudden  changes  in  the  input 
signal  characteristics  due  to  such  events  as 
rocket  ignition,  burnout  and  stage  separation. 
Reference  9  showed  that  the  first  deficiency  could  be 
overcome  by  using  an  "expanding  memory"  filter  algorithm  and 
that  the  second  could  be  overcome  by  applying  multiple  exponen¬ 
tial  filters  having  different  weights  to  the  input  data  and 
using  an  intercomparison  scheme  to  select  the  optimal  outputs. 

The  recommended  second  degree  polynominal  expanding  memory  filter 
is  recursively  updated  as  follows: 

Starting  with  X  =  X°,  X  =  X  =  j  =  0, 

D  =  (j+3)(j+2}(j+l) 

X  =  X+Xt+x't  2/2 

E  =  X°-X  (4.1.2) 

X  =  x+3(3j2+3j+2)E/D 

X  =  X  +  x't  +  18(2j  +  1)E/(Dt  ) 

X*  =  x’  +  60E/(D  t2) 

3  =  3  +  1 

where 

X°  is  the  current  observation, 
t  is  the  time  between  samples, 
j  is  the  sample  number  (j =0,1  ...n)  and 
X,X,  X  is  the  current  state  vector  estimate 
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4.2 


Recommended  Algorithm 
4.2.1  A 1 qor i thm  Development 

The  algorithms  defined  by  Equations  4.1.1  and  4.1.2 
are  capable  of  smoothing  and  differentiating  radar  data  to  form 
an  initial  position/velocity  state  vector  which  can  be  used 
for  generating  real-time  instantaneous  impact  predictions  ( 1 1 P  s )  . 
However,  in  addition  to  smoothing  and  differentiating,  a  range 
safety  digital  filter  must  also  perform  other  important  tasks, 
i ncl udi ng  : 

1.  It  must  edit  the  input  data  and  replace  limited 
quantities  of  detected  bad  data  with  extrapolated 
data. 

2.  It  must  detect  and  adapt  to  sudden  changes  in  target 
motion  or  input  data  quality  by  optimally  modifying 
the  weights  assigned  the  data. 

3.  It  must  estimate  the  uncertainty  (variance)  in  its 
smoothed  position  and  velocity  estimates.  These 
uncertainties  can  be  used  for  generating  IIP 
uncertainty  ellipses  and  for  automatic  best  data 
source  selection. 

The  exponential  filter  algorithm  developed  during  this 
study  is  a  modified  version  of  the  one  developed  in  Reference  9. 
The  previously  developed  filter  had  the  following  features: 

1.  It  accepted  data  from  up  to  four  radars. 

2.  Filter  weights  could  be  selected  to  give  zero, 
first  and  second  degree  polynomial  filtering. 
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3.  "On-Track"  bit  was  utilized. 

4.  Three  sets  of  filter  weights  were  applied  to  each 
channel  of  input. 

5.  A  zero  degree  polynomial  filter  was  used  to  estimate 
noise  variance. 

6.  Adaptation  was  accomplished  through  interfilter 
velocity  comparisons,  through  residual  trend  detection 
or  by  both  methods. 

7.  Limits  could  be  placed  upon  maximum  allowable 
residual  size  and  upon  maximum  number  of  allowable 
extrapolations. 

8.  Initialization  and  reinitialization  of  the  exponential 
filter  was  accomplished  with  the  expanding  memory 
least  squares  polynomial  algorithms  (Equation 

4.1.2  above)  . 

9.  Switching  from  the  expanding  memory  to  the  exponential 
algorithms  was  accomplished  through  comparison  of 
variance  reduction  factors  (VRFs). 

10.  Computed  variances  could  be  modified  with  a  priori 
estimates  and  with  computed  VRFs. 

Most  of  the  above  features  have  been  retained  in  the 
filter  recommended  for  the  P2457  update.  The  most  significant 
modification  recommended  is  that  the  required  subroutine  should 
accept  only  one  channel  at  a  time.  This  will  allow  the  sub¬ 
routine  to  be  more  general  and  will  require  that  calibration 
corrections,  refraction  corrections  and  coordinate  transformation 
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take  place  in  another  subroutine.  The  primary  advantage  of 
this  recommendation  is  that  the  subroutine  will  be  able  to  accept, 
without  modification,  data  from  systems  other  than  land  based, 
monopulse  radars  (e.g.,  telemetered  navigation  and  altimeter 
system  data,  MTRACS  data).  A  flow  chart  of  a  test  bed  version 
of  the  subroutine,  (program  RADAR),  designed  for  evaluating  the 
effectiveness  of  the  various  options,  is  presented  in  Figure 
4. 2. 1.1.  A  listing  of  RADAR,  a  description  of  input  requirements 
and  examples  of  tabular  outputs  are  provided  in  Appendix  C. 

A  description  of  the  major  features  of  the  algorithm  follows. 

4.2.2  Filter  Effective  M emo ry 

The  Reference  9  exponential  filter  weights  (a,6,y)  were 
required  inputs  to  the  program.  In  order  to  make  the  weighting 
more  meaningful  to  users  of  the  subroutine,  the  weighting  in 
RADAR  is  controlled  by  inputting  the  "effective  memory"  of  each 
filter.  The  effective  memory  of  an  exponential  filter  is  here 
defined  to  be  the  number  of  points  in  the  span  of  a  moving  arc 
filter  which  gives  the  same  variance  reduction  in  filtered 
position  data  as  the  exponential  filter  of  the  same  degree,  a 
standard  proposed  in  Reference  10.  The  effective  memory  is  used 
to  compute  the  filter  weights  as  follows: 

Given  the  input  filter  effective  memory,  F  points, 

9  =  1 -4 . 36/ F 
a  =  1-e3 

(4. 2. 2.1) 

8  =  15 ( 1  +  6 ) ( 1-e ) 2 
y  =  lOO(l-e)3 
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The  above  relationship  between  the  movinq  arc  and 
exponential  filter  memories  is  highly  accurate  for  long  memories 
but  only  approximate  for  short  memories. 

4.2.3  Initialization,  Rein i t i a 1 i  z  at i on 

The  recommended  algorithm  for  initializing  the  exponential 
filter  and  reinitializing  after  dropouts  was  given  by  Equation 
4.1.2,  from  Reference  15.  The  initialization  algorithm  is  allowed 
to  expand  its  memory  until  its  position  variance  reduction  factor 
(VRF)  is  less  or  equal  to  the  position  VRF  of  the  desired  expon¬ 
ential  filter.  The  filtering  is  then  transferred  to  the  exponential 
algorithm.  The  VRFs  are  computed  as  follows  (from  Reference  10): 

Position  VRF  (Exponential)  =  (l-e)(l9+248+16S2+6e3 

+  r  4  )/('!  +  1  )  5 

Velocity  VRF  (Exponential)  =  50(l-e)3(49  +  50o  +  1362)/(e-  +  l)5 
Position  VRF  (Expanding)  =  (9Nz+27N+24)/(N+l)3 
Velocity  VRF  (Expanding)  =  ( 1 9200N 2 +74400N+68400 ) / ( N+3 ) 5 

(4. 2. 3.1) 

where  N  is  the  number  of  points  in  the  expanding  memory  filter 
span  and  e  was  defined  in  Equations  4. 2. 2.1.. 

4.2.4  Variance  Estimation 

The  recommended  noise  variance  estimation  scheme  used 
in  RADAR  applies  a  zero  degree  exponential  filter  to  the  squares  of 
the  filter  residuals.  To  avoid  the  use  of  extremely  wild  points 
(which  are  detected  in  the  editing  process),  a  limit  is  set  upon 
the  maximum  residual  allowed  to  be  used.  The  variance  estimation 
algorithms  is  as  follows: 
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V( residual )  =  X°-X 

If  ( ] V  i  ' V  )  ,  V  =  V 
v  1  1  max  max 

S(Variance)  =  S  +  n ( V  2  - S  ) 


(4.2.4. 1) 


where  X°  and  X  were  defined  in  Section  4.1. 

The  variance  estimate  output  by  the  filter  can  be 
increased  by  an  input  a  priori  value  to  account  for  undetectable 
low  frequency  noise.  The  program  multiplies  the  noise  variance 
estimates  by  the  previously  defined  VRFs  (Equation  4. 2. 3.1)  to 
obtain  estimates  of  variance  in  the  smoothed  data. 

4.2.5  Data  Editing 

The  recommended  data  editing  scheme  used  in  RADAR 
compares  the  input  observations  against  observations  computed 
from  an  extrapolated  state  vector.  If  the  differences  (REPS) 
exceed  a  prescribed  value  or  the  offtrack  bit  is  on,  the  input 
observation  is  replaced  with  the  extrapolated  value.  The 
algorithm  proceeds  as  follows: 


c3  ~  C1 


c4  =  3S 


1/2 


If  e4  >ei»  e3 


If  REPS  >^3,  extrapolated  value  is  used 


Here . 


S  was  defined  in  Section  4.2.4  and  t ^  is  an  input 
tolerance  value. 

An  input  value,  NMAX,  places  a  limit  on  the  maximum 
number  of  extrapolations  allowed  before  reinitialization  of  the 
filter  must  take  place. 
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4.2.6 


Filter  Adapyitiqn 

Closely  related  to  data  editing  is  filter  adaptation. 


This  consists  of  varying  the  filter  weights  as  a  function  of 
target  dynamics  or  of  data  quality.  Servo  error  voltage  and 
AGC  measurements  have  been  used  for  this  purpose  in  the  develop¬ 
ment  of  real-time  digital  filters  for  ON-AXIS  radars  (Reference 
15).  In  typical  usage,  if  the  filtered  error  voltages  exceed 
prescribed  values,  the  filters  are  "opened  up",  giving  more 
weight  to  newly  acquired  data  and  reducing  the  effective  filter 
memory.  AGC  measurements  are  used  to  place  limits  on  the  filter 
weights.  Unfortunately,  servo  error  voltage  and  AGC  measure¬ 
ments  are  not  available  at  the  AD  real-time  computers;  alternative 
adaptation  procedures  are  required. 

Two  adaptation  tests  were  developed  during  the  Reference 
9  study.  These  have  been  modified  and  extended  during  this  study. 
The  basic  purposes  of  the  tests  are  to: 

a)  avoid  settling  periods  after  events  such  as  missile 
stage  separation  and  data  dropouts,  and 

b)  detect  possible  filter  bias  or  divergence. 

The  revised  adaptation  process  developed  during  this 
study  and  included  in  the  RADAR  filter  program  makes  use  of  the 
fact  that  several  exponential  filters  can  be  applied  in  less 
computer  time  than  is  required  for  other  types  of  filters. 

This  allows  the  optimum  filter  to  be  selected  through  inter-filter 
comparisons  and  simulates  the  type  of  adaptation  being  accomplished 
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with  error  voltage  and  AGC  measurements  by  other  ranges.  At 
the  same  time,  it  eliminates  response  and  settling  problems 
commonly  experienced  with  other  filters.  In  the  adaptation 
algorithm  developed  for  the  Reference  9  study  and  extended 
during  this  study,  three  exponential  filters  were  applied  to 
each  channel  of  input. 

The  complete  adaptation  process  is  shown  in  the  Figure 
4. 2. 1.1  flow  chart,  beginning  at  entry  point  DD.  Several 
sequential  tests  must  be  passed  by  the  filter  selected  as  best. 

The  tests  include  the  following: 

1)  Tests  to  determine  if  the  filter  is  extrapolating 
due  to  the  editing  procedure  (Section  4.2.5). 

The  longest  filter  which  is  not  extrapolating  is 
normally  designated  best. 

2)  Tests  which  compare  differences  between  filter  outputs 
against  a  tolerance  value  and  count  the  number  of 
successive  differences  which  exceed  the  tolerance. 

A  large  number  of  successive  differences  exceeding 
the  tolerance  is  assumed  to  indicate  that  the  longer 
filter  is  biased  or  diverging,  and  the  shorter  filter 
is  then  designated  best.  The  tolerance  (T0L2) 
and  the  allowed  maximum  number  of  differences  (MAXNV) 
are  input  values. 

3)  Tests  which  count  the  number  of  successive  residuals 
having  the  same  sign.  If  the  number  exceeds  an  input 
value  (NRMAX),  the  filter  output  is  assumed  to  be 
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biased  or  diverging,  and  a  shorter  filter  is 
selected  as  best. 

4.2.7  Filter  C oordinate  Systems 

Reference  9  details  several  advantages  for  filtering 
in  a  system's  natural  coordinate  system  (e.g.,  azimuth, 
elevation,  range)  rather  than  in  transformed  coordinates 
(e.g.,  rectangular,  cartesian).  The  advantages  include: 

a)  Unlike  cartesian  coordinates,  which  are  subject 

to  cross-coupling  of  errors,  the  natural  coordinates 
can  usually  be  considered  independent.  This  simpli¬ 
fies  error  detection  and  analysis. 

b)  Filter  weights  can  be  selected  in  a  more  optimal 
manner.  For  example,  the  FPS-16  range  channel 
normally  requires  a  much  shorter  filter  memory 
than  the  angle  channels. 

The  Reference  9  study  only  addressed  the  use  of  pulse 
radar  data  for  range  safety  control.  Since  that  time,  the  use 
of  additional  systems  has  been  recommended  (e.g.,  telemetered 
navigation  and  altimeter  data,  MTRACS  data).  A  re-examination 
of  the  filter  coordinate  system  problem  during  this  study  has 
led  to  the  conclusion  that  for  the  P2457  update,  it  is  now 
preferable  to  continue  performing  the  digital  filtering  in 
the  test  rectangular  cartesian  coordinate  system.  This  will 
greatly  simplify  the  combining  of  the  pulse  radar  data  with 
data  from  other  systems  (see  Section  5). 
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4.2.8  Output  Character! s  t  i  cs 

Examples  of  the  output  characteristics  of  the  proposed 
adaptive  exponential  filter  are  shown  in  Figures  4. 2. 8.1  and 
4. 2. 8. 2.  The  input  data  was  Site  A-20  radar  range  data  derived 
from  a  simulated  AMRAAM  trajectory  and  realistically  contaminated 
with  noise.  Figure  4.2.8. 1  contain^  nipts  of  the  filter 
velocity  outputs  in  feet/second  versus  time  in  seconds. 

Figure  4. 2. 8. 2  contains  plots  of  the  filter  residual  outputs 
in  feet  versus  time.  The  effective  memories  of  the  shortest, 
intermediate  and  longest  filters  were  21,  41  and  81  points, 
respectively.  The  plots  of  the  outputs  which  the  adaptive 
process  determined  best  at  each  time  point  are  in  the  lower 
right  hand  corner  of  each  figure.  Some  of  the  characteristics 
of  interest  in  the  two  figures  were  the  following: 

1.  Because  the  missile  attained  maximum  acceleration 
within  a  fraction  of  a  second  after  ignition, 
causing  high  order  derivatives  to  appear  in  the  signal, 
all  three  of  the  filters  lagged  to  such  an  extent 

that  the  adaptive  logic  required  them  to  reinitialize 
during  the  first  two  seconds  of  flight  (event  marked 
A  in  the  figures). 

2.  The  41  and  81  point  filters  diverged  at  points 
where  rapid  changes  in  missile  dynamics  occurred. 

The  divergence  was  detected  by  the  adaptive  process 
through  the  trended  residual  test  (see  Section  4.2.6). 
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The  divergence  can  easily  be  seen  in  the  residual 
plots  of  Figure  4. 2. 8. 2  at  the  points  marked  B. 

The  effect  on  the  velocity  outputs  can  also  be  seen 
at  the  points  marked  B  in  Figure  4. 2. 8.1. 

3.  The  effectiveness  of  the  multiple  filter  adaptive 
process  can  be  seen  in  the  "Best  Filter"  output 
plot  of  Figure  4. 2. 8.1.  Note  that  the  adaptive 
process  tended  to  select  the  shortest  memory  filter 
output  during  periods  of  high  dynamics  and  the 
smooth  longer  memory  filter  outputs  during  periods 
of  low  dynamics. 
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Residual  Output 

21  Point  Effective  Memory  Filter 


Residual  Output 

41  Point  Effective  Memory  Filter 
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REAL-TIME  TELEMETRY  DATA  PROCESSING  ALGORITHMS 


5.1  Background 

Many  of  the  advanced  weapon  programs  planned  for 
development  by  AD  will  require  low  level  tests  in  which  the 
acquirpr!  radar  data  will  be  severely  degraded  due  to  clutter, 
multipath  and  refraction  effects.  Because  of  the  high  dynamics 
of  the  vehicles,  radar  track  will  be  difficult  to  reacquire 
if  it  is  lost.  An  alternative  to  the  use  of  only  radar  data 
for  range  safety  control  of  advanced  weapon  tests  is  the  use 
of  telemetered  on-board  sensor  data  as  back-up  to  the  radar 
data  or  in  a  combined  radar/telemetry  (TM)  data  trajectory 
solution.  On-board  sensor  data  which  might  be  used  to  augment 
the  tracking  radar  data  include: 

•  Inertial  navigation  system  (INS)  data 

•  Radio  navigation  system  (RNS)  data 

•  Radar  or  barometric  altimeter  data 

•  Status  information  (e.g.,  seeker  lock-on  indication, 
loss  of  command/control  signal) 

One  of  the  Task  SER  82-1  conclusions  was  that: 
"Algorithms  for  the  use  of  telemetered  inertial  navigation  and 
altimeter  data  should  be  included  in  the  P2457  update".  This 
conclusion  was  based  largely  upon  the  analyses  conducted  during 
Task  SER  81-3  (Reference  11).  One  of  the  conclusions  resulting 
from  81-3  was  the  following: 

Vertical  coordinate  and  IIP  data  can  be  improved  in 

quality  by  an  order  of  magnitude  (and  sometimes  more) 
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by  judiciously  combining  the  radar  data  with  tele¬ 
metered  altimeter  data  and  INS  position  and  velocity 
data.  The  combinative  methods  developed  during  this 
study  are  both  feasible  and  desirable  for  the  real¬ 
time  estimation  of  present  position  coordinates  and 
I  IPs  for  use  in  the  range  safety  control  of  AD  weapon 
tests.  The  combinative  methods,  involving  the 
application  of  least  squares  exponential  filters  to 
radar-INS  data  differences  and  to  radar-altimeter 
data  weighted  averages,  are  described  in  this  report 
(Reference  11). 

An  example  of  the  improvement  in  IIP  quality  which  can 
be  achieved  by  combining  radar  and  TM  data  is  shown  in  Figures 
5.1.1  through  5.1.4,  the  results  of  a  realistic  simulation  of 
data  acquired  on  a  low-level  test  over  the  Gulf.  Figure  5.1.1 
shows  the  noisy,  but  unbiased,  IIP  error  trace  resulting  from 
the  use  of  radar  data  only.  Figure  5.1.2  exhibits  a  typical 
smooth,  but  trended,  IIP  error  trace  resulting  from  INS  data 
only.  Figure  5.1.3  shows  that,  by  combining  the  radar  and  INS 
data,  a  relatively  smooth,  untrended  IIP  trace  results. 
Similarly,  Figure  5.1.4  shows  that  by  combining  low  elevation 
radar  data  with  altimeter  data,  the  IIP  trace  can  be  greatly 
improved . 

A  purpose  of  this  study  was  to  examine  existing  radar/ TM 
combinative  algorithms,  develop  new  ones  if  necessary,  and 
document  the  recommended  algorithms  in  a  form  conducive  to 
efficient  incorporation  in  a  P2457  update. 
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5.2  Coordinate  Systems 

There  is  an  enormously  wide  variation  in  the  types  of 
coordinate  systems  in  which  telemetered  inertial  navigation 
and  guidance  system  data  have  been  presented.  Variations 
include  the  following  types  of  coordinates: 

I  Geodetic  latitude,  longitude,  altitude 
t  Rectangular  cartesian 
•  Scrunched  ( non-orthogonal  axes) 

0  Eartti  fixed 
I  Inertial 
0  Topocentric 
0  Earth  centered 
0  Flight  line  oriented  axes 
0  North-south,  east-west  oriented  axes 
0  Tilted  axes 

A  common  additional  problem  in  such  data  is  that  the 
position  data  may  be  degraded  by  telemetering  an  inadequate 
number  of  bits  to  provide  a  smooth  trajectory  plot  (e.g.,  the 
least  bit  telemetered  by  one  aircraft  INS  is  .001°  in  latitude 
and  longitude,  or  as  much  as  365  ft.)  This  deficiency  can 
often  be  overcome  by  numerically  integrating  the  telemetered 
velocity  data,  for  which  the  least  bit  is  normally  less  than 
one  ft/sec. 

An  in-depth  analysis  of  all  of  the  coordinate  system, 
least  bit  and  other  problems  involved  in  transforming  telemetered 
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INS  data  into  a  form  suitable  for  use  in  range  safety  control 
was  considered  to  be  beyond  the  scope  of  this  study.  These 
problems  must  be  considered  on  an  individual  weapon  system 
basis  and  a  suitable  transformation  subroutine  developed  when 
a  requirement  for  real-time  use  in  range  safety  control  is 
levied.  It  will  be  assumed  here  that  the  telemetered  data  has 
been  transformed  into  a  form  suitable  for  input  to  the  XYZLLH 
subrout i  ne. 

5.3  Recommended  R_adar/INS  Combinative  Algorithm! 

It  was  concluded  in  Task  SER  81-3  (Reference  11) 
that  a  Kalman  or  batch  filter  "best  estimate  of  trajectory 
(BET)"  type  approach  to  combining  radar  and  TM  data  in  real¬ 
time  was  too  time  consuming  for  use  on  the  multiple  vehicle 
test  scenarios  planned  for  the  Eglin  ranges.  A  re-examination 
of  the  problem  during  this  study  has  reaffirmed  this  conclusion. 

It  is  recommended  that  the  algorithms  developed  in  Task  81-3, 
modified  to  utilize  the  filter  developed  during  this  study,  be 
included  in  the  P2457  update.  The  complete  algorithm  can  be 
expressed  as  follows: 

X  B  X INS  +  ^XR  "X I N S ^  (5.3.1) 

where 

X  R  =  Best  estimate  of  the  position/velocity/acceleration 
vector,  using  radar/INS  combination 

X  i n $  =  Vector  of  inertial  navigation  indicated  coordinates 

X„  =  Vector  of  radar- i nd i cated  coordinates 
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|  Here  the  bar  represents  low  pass  filtering  of  the 

radar-INS  differences.  It  is  recommended  that  the  adaptive 
1  exponential  filter  developed  in  Section  4  of  this  report  be 

used  to  perform  the  required  filtering.  Because  the  filtering 
process  of  Equation  5.3.1  is,  in  essence,  a  method  for  estimating 
low  frequency  errors  in  the  INS  data  (e.g.,  the  quadratic  error 
trend  of  Figure  5.1.2),  the  required  effective  memory  inputs 
to  the  filter  will  be  much  longer  than  for  radar  data  alone. 

A  minimum  effective  memory  of  10  seconds  is  suggested. 

A  conceptual  flow  chart  of  the  combinative  radar/INS 
data  processing  is  shown  in  Figure  5.3.1. 

5.4  R eco mm e n ded  Radar/Al ti m ete r  Co m b i n a t i v e  A 1 g or i t hm 

Since  radar  or  barometric  altimeters  supply  only  one 
of  the  three  coordinates  required  to  define  a  present  position 
vector,  a  different  combinative  algorithm  than  that  used  in 
combining  radar  and  INS  data  is  required.  The  combinative 
algorithm  recommended  here  is  a  modified  version  of  the 
algorithm  proposed  in  Reference  11,  the  only  major  modification 
being  in  the  estimation  of  the  weights  to  be  used  in  the 
weighted  average  process.  The  weighted  average  algorithm 
used  to  obtain  a  best  estimate  of  the  Y  (up)  coordinates  is 
as  follows: 

Yb=  (°y2  +  °A  ^  ^  (YR  °y2  +  Ya  °A2)  (5.4.1) 

R  R 
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Conceptual  Flow  of  Radar/INS  Combinative  Data  Process 


Figure  5.3.1 
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where 


Y  R  =  weighted  least  squares  estimate  of  Up  coordinate 
°  relative  to  the  radar 

Yd  =  radar  indicated  Y  ,  corrected  for  refraction 
K  2  2 

Y^  ~  Ha  -  Hr  -  R  cos  E/2Re,  a  good  approximation 
of  altimeter  indicated  Y  for  AD  weapon  tests 

H  R  =  altitude  of  radar  site 

=  altimeter  indicated  altitude  above  sea  level 

R  =  mean  radius  of  earth 
e 

2  2 

Rather  than  using  a  priori  values  for  ov  and  o  ,  as 

yR  A 

was  suggested  in  Reference  11,  it  is  recommended  that  these 
variances  be  estimated  by  passing  Y^  and  through  the  exponential 
filter  and  obtaining  variance  estimates  by  Equations  4. 2. 4.1. 

After  obtaining  the  new  estimate  of  Yg  using  Equation  5.4.1, 
improved  estimates  of  X  and  Z  can  be  computed  as  follows: 

XB=  (r2  .  Y2)1/2  Cos  A  (5.4.2) 

Zg  =  (R2  -  Y2)1/2  Sin  A 

where  A  is  the  radar  azimuth  measurement  and  R  is  the  refraction 
corrected  radar  range. 

After  obtaining  the  improved  estimates  of  X,  Y,  Z  using 
the  combinative  algorithm,  the  data  should  be  input  to  the 
exponential  filter  to  obtain  a  smoothed  state  vector  and  the 
related  uncertainty  estimates.  A  conceptual  flow  chart  of  the 
radar/altimeter  data  combinative  process  is  shown  in  Figure 
5.4.1. 
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Conceptual  Flow  of  Combinative  Radar/Altimeter  Data  Processing 

Figure  5.4.1 
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APPENDIX  A 


LAGRAN  Interpolation/Extrapolation  Subroutines 

Subroutine  LAGRAN  performs  three  variable, 

2nd  degree  polynominal,  Lagrangian  interpolation  of  tables 
of  the  form  shown  in  Figure  A.l.  A  FORTRAN  listing  of  the 
subroutine  is  presented  on  Paqes  A-3  through  A-5.  The 
elements  of  the  calling  sequence  are  as  follows: 

XX  is  the  X  coordinate  of  the  desired  interpolated 
or  extrapolated  point, 

YV  is  the  Y  coordinate  of  the  desired  interpolated 
or  extrapolated  point, 

11  is  the  Z  coordinate  of  the  desired  interpolated 
or  extrapolated  point, 

X  is  the  array  of  X  coordinates, 

Y  is  the  array  of  Y  coordinates, 

Z  is  the  array  of  Z  coordinates, 

TAB  is  the  three  dimensional  array  of  function 

values,  given  by  TAB  (I,J,K)  =  F(X(I),  Y(J),Z(K)), 

NX  is  the  integer  number  of  elements  in  the  X  array, 

NY  is  the  integer  number  of  elements  in  the  Y  array, 

NZ  is  the  integer  number  of  elements  in  the  Z  array, 

M  is  3  for  second  degree  polynominal  interpolation, 

IEXTX  is  an  integer  returned  by  the  subroutine 
which  has  the  following  meanings: 

IEXTX  =  1  if  extrapolation  occurred  above  X  table 
IEXTX  =  0  if  interpolation  occurred 
IEXTX  =  -1  if  extrapolation  occurred  below  X  table 
IEXTX  =  M  if  M>N Y 

IEXTY  is  the  same  for  the  Y  array  as  IEXTX, 

ANS  is  the  returned  interpolated  or  extrapolated 
function  approximation. 
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LAC, RAN  Subroutine  HOkTRAN  LISTING 


SUBROUTINE  LAGRAN(XX,YY,ZZ,X,Y,Z.TAB,NX,NY,NZ,M,1EXTX,IEXTY,ANS) 
C 

C  THREE  VARIABLE  L AGRARIAN  INTERPOLATION 
C 

DIMENSION  TAB(NX,NY,NZ),X(n,Y(l).Zm,XLAG(100).DLAG(3) 

LOGICAL  QUITX.GUITY 
N=NZ  -  1 
DO  29  LL=2.N 

IF ( ZZ  -  Z(LL+1 ) >30.29,29 

29  CONTINUE 
LL=NZ  -  1 

C 

C  INITIALIZE 
C 

30  GUI TX=. FALSE. 

GUI TYr. FALSE. 

Ml  =M 

I£XTX=0 

IEXTY=0 

C 

C  FIND  THE  RANGE  QF  INTERPOLATION  ALONG  X 
C 

IF ( Ml  .LE.  NX) GO  TO  1 

Ml -NX 

IEXTX=M1 

1  DO  2  1=1. NX 

IF(XX  -  X(I) )3.9.2 

2  CONTINUE 
IEXTX=1 

GO  TO  8 

3  MIDX=I 

IF(MIDX  .GT.  1)G0  TO  4 

IEXTX=-1 

GO  TO  5 

4  IF(ABS(XX  -  X(MIDX) )  .GE.  ABS(XX  -  X(MIDX  -  1)))G0  TO  5 
ISXPT=MIDX  -  (Ml  /  2) 

GO  TO  6 

5  ISXPT=MIDX  -  ((Ml  +  1)  /  2) 

G  IFdSKPT  .GT.  0)G0  TO  7 

ISXPT=1 

7  IEXPT=ISXPT  ♦  Ml  -  1 

IFdEXPT  .LE.  NX) GO  TO  10 

8  ISXPT=NX  -  Ml  +  1 

IEXPT=NX 

BO  TO  10 

9  GUITX=.TRUE. 

C 
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LAGKAN  Subroutine  FORTRAN'  LISTING 


C  FIND  THE  RANGE  OF  INTERPOLATION  ALONG  Y 
C 

10  IF (Ml  .LE.  NY) GO  TO  11 
M1=NY 

IEXTY=M1 

11  DO  12  JU  .NY 

IF ( YY  -  YU)  >13,19,12 

12  CONTINUE 
I  EXT Y= 1 

GO  TO  18 

13  MIDY=J 

IFIMIDY  .GT.  1 ) GO  TO  14 

IEXTY=-1 

GO  TO  15 

14  IF ( ABS( YY  -  Y(MIDY) )  .GE.  ABS( YY  -  Y01IDY  -  1)))G0  TO  15 
ISYPT ~MIDY  -  (Ml  /  2) 

GO  TO  16 

15  ISYPT=MIDY  -  ((Ml  +  1)  /  2) 

16  IFdSYPT  .GT.  0)G0  TO  17 
ISYPT=1 

17  IEYPTUSYPT  +  Ml  -  1 

IF ( IEYPT  .LE.  NY) GO  TO  20 

18  ISYPT=NY  -  Ml  +  1 
IEYPT=NY 

GO  TO  20 

19  GUITY= .TRUE. 

20  N=0 

DO  28  LLL=(LL  -  1),<LL  +  1) 

N=N  +  1 

IF(GUITX  .AND.  QUITYJGO  TO  27 
IF(QUITX  .AND.  .NOT.  QUITY)GO  TO  24 
IF( .NOT.  GUITX  .AND.  3UITY)G0  TO  23 
C 

C  INTERPOLATE  ALONG  X 
C 

21  DO  22  K=ISYPT, IEYPT 

22  CALL  INTER(XXrTAB( 1 rK.LLL) .X. ISXPT , IEXPT.XLAG(K) ) 

GO  TO  26 

23  CALL  INTER(XXf TABU r JfLLL) ,X, ISXPT » IEXPTrXLAQ( J) ) 
DLAG(N)=XLAG( J) 

GO  TO  28 
C 
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LAGRAN  Subroutine  FORTRAN  LISTING 


C  INTERPOLATE  ALONG  Y 

C 

24  DO  25  L= ISYPT. IEYPT 

25  XLAG(L) rTAB( I >L,LLL) 

2G  CALL  INTER(YY,XLAG.Y, ISYPT, IEYPT, YLAG) 

DLAG(N) - YLAG 
GO  TO  28 

27  DLAG(N) =TA3(I,J,LLL) 

28  CONTINUE 

C 

C  INTERPOLATE  ALONG  Z 

C 

B 1 =  Z ( L L  -  1)  -  Z(LL) 

B2=Z(LL  -  1)  -  Z(LL  +  1) 

B3=Z (LL)  -  Z(LL  +  1) 

A1  =ZZ  -  Z ( LL  -  1) 

A2=ZZ  -  Z(LL) 

A3=ZZ  -  Z(LL  +  I) 

Cl =DLAG < 1 )  /  ( B 1  *  B2) 

C2=-DLAG(2)  /  (BI  »  B3) 

C3=DLAG(3)  /  (B2  *  B3) 

ANS=A2  *  A3  *  Cl  +  A1  *  A3  *  C2  +  A1  #  A2  *  C3 

RETURN 

END 


SUBROUTINE  INTER(PT,FUN,ORD, ISPT, IEPT.ALAG) 
C 

DIMENSION  FUN(l) -QRD(l) 

SUM=0. 

PROD=l . 

DO  3  J=ISPT , IEPT 

PROD=PROD  *  (FT  -  ORD(J)) 

DENOM=l . 

DO  2  K=ISPT , IEPT 
IF(K  .NE.  J)GO  TO  1 
D=PT 
GO  TO  2 

1  D=ORD<  J) 

2  DENOM=DENOM  *  (D  -  ORD(K)) 

3  SUM=SUM  +  FUN(J)  /  DENOM 
ALAGrSUM  *  PROD 

RETURN 

END 
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1 


-  -WS . *  COS  ( WD . ) 

1  i 

v;  r. .  -  -  w  s  .  *  sin  (  w  d  . ) 

ii  i 

w s  .  :  wind  speed  at  altitude  h . 

i  i 

w a  =  wind  direction  at  altitude  h. 

i  i 

The  c  nr  pc  nents  'X..  '7.  are  sunmed  over  all  intervals 

i  i 

b  e  t  w  (;  e  n  the  ground  and  the  altitude  at  w  hi  ch  t  h  e  parachute  o  p  e  n  s 
to  give  the  total  northward  and  eastward  drifts  of  the-  parachute 
at  impact.  Subroutine  XYZLLH  can  then  be  entered  to  determine 
the  i  it:  pact  point  latitude  and  longitudes.  A  FORTRAN  lasting 
of  a  version  of  the  above  parachute  impact  predictor,  which  is 
included  in  the  DENG1  program  (Reference  14),  is  presented  in 
n  a  g e s  B  -  3  and  B  -  4  .  The  linear  i  n  t e rpol  a  ti  on  s u b r t . u  t  i  n e  tall-,  d 
by  the  predictor  is  listed  on  Page  B-5. 
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r  ACT  i'EHLICTOR  (PIP) 


SUBROUTINE  PIP( ICOL, JP IP) 

C—  PARACHUTE  IMPACT  PREDICTOR - 

P < 2 7 1 COL )  IS  X  NORTH, FT. 

P( 3, ICOL)  IS  Y  UP,  FT. 

P ( 4 , 1  COL )  IS  Z  EAST,  FT. 

C  Pill, ICOL)  IS  ALTITUDE,  FT. 

REAL  LAMQRI 

COMMON  /TABLES/  WEATHN ( 50 , 5 ) , WEATHER (5,50), AREAN , I  WE ATH , I  CD , P ( 1 1 , 7 
+  >  > El ,E2, IFLAG, DT , PHI , RLAMDA , ANS ( 10) , PLAT , PLCN, PSI , FL ,H,CFL,SFL 
COMMON  /AERO/  Q , CMACH , TMACHt 1 5 ) , TAL PHA ( 1 0) , TA3CNA (100), TABCA (100), 
1TABTHR (30) , NALPH , NP , CN , CA , NMACH , AREFS , CAERO, TM ACH1 (15), TABCL (100), 
2TCA( 100) ,NMACH1 ,NCA, CDS'/ST (40) ,NCDP,SPINAT 
COMMON  /MREFD/  AREF ,  WE  I GHT , AMASS , CDO , CLO , CDA , CLA . LQD , ABET A , GLATLM , 
+  WQQT ( 10) ,ALOSS 

COMMON  /TRAJ/  PUHIST  <  7 , 50 ) ,PV( 7) , POT ( 7 ) , PH  I  OR I ,LAMORI , PS  I  OR I 
COMMON  /FLAGS/  NOWIND, NOLIFT, NODRAG, NOTHRST, IOUT , LCOUNT, SPINFLG, 

+  J IMMY, FACTOR , I  ENTRY , I  CHUTE , GODOWN, II4ING  > WINSPFS 
COMMON  /BDFG/  BADFLG,EDALT,HTLIM, ALTSTP, TIPIP, TFLMAX 
DIMENSION  WN(50) ,WE<50> ,0T<50> , ALT(50) 

FETID  =  364566.0 

DTR  =  0.0174532925 

NALTS  =  50 

NLCDS  =  NCDP  *  2 

BETA  =  WEIGHT/CDSVST(NLC0S) 

4  DO  5  J=5, 10 
5  P( J, ICOL)  =  0.0 
DO  10  1=1, NALTS 
NWTH  =  1-1 

IF ( WE A i HER (2, I) .EQ.O.O)  GO  TO  12 
WD  =  WEATHER(3, I )*DTR 
WS  =  WEATHER (4,1) 

WN( I )  »  -WS  *  COS(WD)  #  1.S877 
WE( I )  =  -WS  *  SIN(WD)  *  1.6877 
D  =  WEATHER(2, I ) 

UT( I )  =  SORT ( 103075B .0  *  BETA/D) 

10  ALT ( I )  =  WEATHER (1,1) 

12  NALTS  =  NWTH 
Z  =  P ( 1 1 , ICOL ) 

CALL  INPOL(ALT ,WN,Z,WY, NALTS) 

CALL  INPOL(ALT, WE, Z,WX, NALTS) 

CALL  INPOL( ALT , VT ,Z,VTZ , NALTS) 

SUMTIM  =  P( 1 , ICOL) 

XIP  «  0.0 
YIP  =  0.0 
DO  30  I =Z, NALTS 
IF(ALT(I)  -  Z)  15,20,20 
15  DALT  =  ALT ( I )  -  ALT(I-l) 

DELT  =  2.0  *  DALT  /  (VT(I)  +  VT(I-l)) 

XIP  =  XIP  +  0.5  *  DELT  *  (WE(I)  +  WE(I-l)) 
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FORTRAN  LISTING  OF  PARACHUTE 
IMPACT  PREDICTOR  (PIP) 


YIP  =  YIP  +  0.5  *  DELT  *  <WY  +  WN'(I-l)) 

SUMTIM  =  SUMTIM  +  DELT 
GO  TO  40 
30  CONTINUE 

40  IF< SUMTIM  .GT.  0.  .AND.  SUMTIM  .LE.  TFLMAX1G0  TO  45 
BADFLG-1 . 

GO  TO  9000 

45  P(2, I COL)  =  P(2 , 1  COL)  +  YIP 
P(3, I COL)  =  0.0 
P(4 , 1 COL)  =  P  <  4  r I COL)  +  XIP 
P(l.ICOL)  =  SUMTIM 
P(  11 , I COL)  =  0.0 
PHRAD  =  PHIORI  *  DTR 
COPHI  =  COS(PHRAD) 

JPIP  =  1 

PARLAT  =  PHIORI  ♦  (P(2,IC0L)/FETID) 

PARLON  =  L AMOR I  -  (P(4, ICOL)/!FETID*COPHI > > 

WRITE(G,900)  SUMTIM, VT<1) 

900  FORMAT ( IX, 'IMPACT  TIME  (SEC)  AND  VELOCITY  (FT/SEC) ' ,ZX,2(3X,F8. 2) ) 
WRITE <6, 1000)  PARLAT, PARLON 

1000  FORMAT! IX, 'PARACHUTE  GEODETIC  POSITION  AT  IMPACT' ,2<2X,F12. 7) ) 
WRITE(E,1010)  (P( JJ,ICOL) , JJ=2>4) 

1010  FORMAT! IX, 'PARACHUTE  CARTESIAN  POSITION  AT  IMPACT' ,3! 1X.F1Z. 1 ) ) 
PIND  *  1. 

WR I TE ( 1 0 , 4450 )  SUMTIM, PARLAT , PARLON , P I ND 
4450  FORMAT <1X,F12.4,2(2X,F10.7),F4.1) 

9000  RETURN 
END 
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FORTRAN  LISTING  OF  LINEAR 
INTERPOLATION  SUBROUTINE  (INPOL) 


SUBROUTINE  INPOL<X,Y,T,FT,NALTS) 

PT MENS I ON  X(50) .  Y(50) 

IF(X(1)-T)  11.11.4 
11  DO  1  1=1 , NALTS 
I F ( X  < I )-T)  1,2,3 

1  CONTINUE 
GO  TO  4 

2  FT  =  Y(I) 

GO  TO  4 

3  K  =  1-1 

FT  =  Y<X)  +  <Y( I )-Y(K) )  *  (T  -X(K))  /  (X(I)-X(K)) 

4  RETURN 
END 
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Adaptive  Exponential  Radar  Data 
Filter  Program  (RADAR) 

A  FORTRAN  listing  of  a  test  bed  adaptive  exponential 
radar  data  filter  program  (RADAR)  developed  during  this  study 
is  presented  on  Page  C-3  through  C-7.  Required  inputs  to  the 
program  are  as  follows: 

Control  Constants 


Name 

Meaning 

START 

Start  time  (seconds) 

STOP 

Stop  time  (seconds) 

NMAX 

Maximum  number  of  extrapol ati ons  allowed 
before  reinitialization  occurs 

RESMAX 

Maximum  size  of  residual  to  be  used  in 
computing  noise  variance  estimates 

SWATE 

Alpha  filter  weight  to  be  used  in  computing 
variance,  noise  variance 

MAXNV 

Maximum  number  of  successive  interfilter 
velocity  differences  allowed  to  exceed  T0L2 

NRMAX 

Maximum  number  of  successive  residuals  of 
some  sign  allowed 

IGA  (6  values) 

Print  output  every  I GAP  point 

NPOS 

Word  position  on  tape  of  measurement  to  be 
filtered 

TSCALE 

Time  scale  factor 

TOL 1 

Discrepancy  vector  tolerance;  used  to 
determine  need  for  extrapolation 

C-l 


Name 


Meaning 


T0L2  Tolerance  on  differences  between  velocity 

estimates  of  two  filters 

FILLEN  (3  values)  Effective  memory  (points)  in  shortest, 

middle  and  longest  filter 

SPRIOR  Estimated  standard  deviation  of  low 

frequency  noise  component 

NPLOTT  Plots  desired:  YES  or  NO 


Tape  Input:  The  program  RADAR  currently  accepts  only  the 
output  of  the  AMRAAM  radar  data  simulation 
program  SITE. 

Examples  of  the  printed  ouput  of  RADAR  are  presented  in  Figures 
C.l  through  C.6. 
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PROGRAM  RADAR!  INPUT .OUTPUT , TAPE2, TAPE3, TAPE4 . TAPE5,  TAPEG, 

1  TAPE15= INPUT , TAPE  10) 

DIMENSION  XPP ( 4 ) . XVP < 4 ) . XAP ( 4 ) , REPS ( 3 ) , NEX ( 3 ) . I FLAG! 3) , FWATE(3>3) 
DIMENSION  RESID(4)  ,RNSIG(3)  ,RPSIG(4)  .RSSIG(3,2)  ,VRF<3.2)  ,VRFE(4,2) 
DIMENSION  NP( 3 ) . PRES ID( 3 ) .NRESI D( 3) , THETA! 3 ) .IPRINT(G) 

DIMENSION  I GAP ( S) , WORD! 72),FILLEN(3),KFG(12),TIT(3),LFLAG(3) 

CHARACTER  *3  NPLOTT 

REWIND  10 

NV32=0 

NV21 =0 

DO  990  1  =  1,  S 
990  I  PRINT ( I ) =-1001 
DO  1010  JR= 1 ,3 
RESID( JR ) =0. 

PRESID! JR ) =0. 

RNSIG! JR) =0. 

REPS( JR) =0. 

NRESID! JR) =0 
NP( JR) =0 
NEX( JR)=0 
I FLAG ( JR ) = 1 
XPP(JRJ=0. 

XVP! JR) =0. 

XAP( JR)=0. 

DO  1010  K  =  1 ,2 
VRFE! JR,K)=0. 

1010  RSSIG! JR,K)=0. 

READ! 15, *)START, STOP, NMAX,RESMAX,SWATE,MAXNV,NRMAX, IGAP.NPDS, 

1  TSCALE 
C 

C  DISCREPANCY  AND  VELOCITY  TOLERANCE 

C 

READ ( 15 , * ) T0L1 r T0L2 

INPUT  FILTER  LENGTHS 

READ(15,*)(FILLEN!JR),JR=1,3) 

DO  G  JR=1,3 

THETA(JR)=1  -  4.3G  /  FILLEN(JR) 

FWATE(JR»1)=1  -  THETA(JR)*#3 

FWATEl JR,2)=15  *  <1  +  THETA(JR))  *  (1  -  THETA( JR) )**Z 
FWATEl JRr3)=100  *  (1  -  THETA! JR) )*«3 

PRIOR  ESTIMATE  OF  OBSERVATION  NOISE 

READ! 15f *)SPRIOR 
SPRI02=SPRI0R**2 

INPUT  PLOT  DETERMINATION 
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PEAD< 15.  * )NPLOTT 

PRINT  11,START,ST0P,NMAX,RESMAX,SWATE,MAXNV,NRMAX,IGAP,NP0S, 

I  T SCALE. T0L1 , T0L2. < <FWATE< IR, JR) , JR=1 ,3) . IR=1 ,3) , SPRIOR, NPLOTT 
11  F0RMAT(1H1, BOX, 'INPUT  CONSTANTS '/// 10X, 'RAW  START  TIME ' , T40.F1 2.2/ 
1/1  OX » 'RAW  STOP  TIME', T40.F12.2//10X,  'MAX  EXTRAPOLATIONS  ALLOWED'. 
2T40. I5//10X. 'MAX  RESIDUAL' .T40.F10.3//10X, 'VARIANCE  FILTER  WEIGHT' 
3, T40.FS.3//10X. 'MAXIMUM  NUMBER  OF  VELOCITIES', T40, I3//10X, 

4 'MAX  RESIDUALS  ALLOWED T40, 15// 

GlOX, 'PRINT  TIME  INTERVAL' .T40.GI5// 

710X, 'WORD  POSITION  TO  BE  FILTERED ', T40, 12// 

91 OX, 'TIME  SCALE  FACTOR' ,T40,F10.4//10X, 'DISCREPANCY  TO!  FRANCE', 
*T40, F10.5//10X, 'VELOCITY  TOLERANCE', 

*T40,F10.5//10X, 'FILTER  WEIGHTS ' ,T40, 

13(3F10.7/T40)//10X, 'OBSERVATION  NOISE  ESTIMATE ', T40, 

2F10.5//10X, 'PLOTTING  DETERMINATION ' ,T40,A3/lHt ) 

DO  1080  JR=1 ,3 
TH=THETA( JR) 

0NE=1.  /  (TH  +  1 . )**5 
RPSIG( JR) -SPRI02 

VRF ( JR,  1 )  =  ( 1 .  -  TH)  *  ONE  *  (19.  +  24.  *  TH  +  18.  *  TH**2  +  8.  * 
1  TH**3  +  TH**4) 

1080  VRF( JR,2)  =  ( 1 .  -  TH)**3  *  50.  *  ONE  *  (49.  +  50.  *  TH  +  13.  * 

1  TH**2) 

4  READ( 10) TIME,  (WORD(N) ,N=1 ,72) , (KFG(N) ,N=1,12),(TIT(N),N=1,3) 

IF ( ABS(TIME  -  START)  .LT.  .001 )G0  TO  7 
GO  TO  4 

7  IF ( (TIME  -  STOP)  .GT.  .001 )G0  TO  999 
KFLAG-KFG< (NPOS  -  1)  /  8  ♦  1) 

RO=WQRD(NPOS) 

BESF=0. 

DO  GO  JR=1 ,3 
T0L3=T0L1 

T0L4=3.  *  SORT (RPSIG( JR) ) 

IF(T0L4  .GT.  T0L1 )TOL3=TOL4 

XPP( JR)=XPP(JR)  +  .1  *  XVP(JR)  +  .005  *  XAP(JR) 

XVP( JR)=XVP( JR)  +  .1  *  XAP(JR) 

IF(KFLAG)905,905, 908 
905  REPS(JR)=0. 

GO  TO  52 

908  IF( IFLAG( JR) )910,51,910 

910  IF(NP( JR)  .EQ.  0)XPP(JR)=R0 

51  REPS (JR) =R0  -  XPP(JR) 

IF(IFLAG( JR) >915,920,915 

915  LFLAG(JR)=1 

IF(NP( JR)  .LT.  5)LFLAG( JR)=0 
IF(NP( JR)  -  10)80,80,920 
920  IF(ABS(REPS( JR) )  -  T0L3)9?5,925,52 

925  NEX( JR)=0 

GO  TO  60 

52  NEX( JR)=NEX( JR)  +  1 
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IF (NEX( JR)  -  NMAX)60,S0,54 
54  1FLAG( JR) =1 

NEX( JR  >  =0 
NP( JR) =0 
GO  TO  910 
SO  CONTINUE 

DO  80  JR= 1 f 3 

IF(IFLAG< JR) >90,70,80 

70  IF(NEX( JR) >80,71 ,80 

71  XPP(JR)=XPP< JR)  +  FUATE(JR,1)  *  REPS(JR) 

XV P ( JR ) =XUP( JR)  +  FWATE( JR, 2)  *  REPS(JR) 

XAP( JR) -XAP( JR)  +  FWATE( JR, 3)  *  REPS(JR) 

80  CONTINUE 

DO  100  JR=1 ,3 

IF ( I FLAG ( JR>  .EG.  0)G0  TO  100 
XNP=NP< JR) 

IF( (XNP  .GT.  10.)  .AND.  <NEX(JR)  .GT.  0))G0  TO  95 
DEIFREPS(JR)  /  ((XNP  +  3.)  *  (XNP  +  2.)  *  (XNP  +1.)) 

XPP( JR) =XPP( JR)  +  (9  *  XNP**2  +  9  *  XNP  +  6)  *  DEN 
XVP( JR)=XUP( JR)  +  (3S0.  *  XNP  +  180.)  *  DEN 
XAP( JR)=XAP( JR)  +  SOOO.  *  DEN 

95  VRFE( JR , 1 ) = (9.  *  XNP**2  +  27.  *  XNP  +  24.)  /  (XNP  +  l.)**3 
VRFE( JR, 2) =(19200.  *  XNP**2  +  74400.  *  XNP  +  GS400.)  / 

1  (XNP  +  3. )**5 

100  CONTINUE 

DO  120  JR-1 ,3 
RESID( JR) =R0  -  XPP(JR) 

TEMP-RESID( JR) 

IF (ABS(TEMP)  -  RESMAX)114, 114,112 
112  TEMP=RESMAX 

114  RNSIG( JR) =RNSIG( JR)  +  SWATE  *  (TEMP**Z  -  RNSIG(JR)) 

RPSIG( JR) =ABS(RNSIG( JR) )  +  SPRIOZ 
120  CONTINUE 

DO  140  JR=1,3 
DO  140  K=1 rZ 

IF ( IFLAG( JR) ) 130,124,130 
124  RS5IG( JR,K)=VRF( JR,K)  *  RPSIG(JR) 

GO  TO  140 

130  RSSIG( JR,K)=VRFE( JRfK)  *  RPSIG(JR) 

140  CONTINUE 

DO  160  JR=1»3 

IF(IFLAG(JR))144, 160,144 
144  IF(VRF(JR,1)  -  VRFE( JR, 1) )148, 152,152 

148  NP( JR)=NP( JR)  +  1 

GO  TO  160 
152  NP( JR)=0 

IFLAG( JR)=0 
160  CONTINUE 

DO  175  JR=1 ,3 

IF (RESID( JR)  *  PRESID(JR) >170, 170,168 
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1G8  NRESID<  JR)=N-T5ID(.JR)  +  1 

GO  TO  175 

170  NREISID<  JR) -0 

175  PRES! D<  JR) -RESIDf  JR) 

IF(LFLAG(3)  >200-  200,176 
17G  I F ( ORF E ( 3 1 1 )  -  0R(-E<2, 1 ))  177, 177.201 

177  IF ( VRFE( 3,1)  -  VRFE(1 , 1 ) >181 , 181 ,Z10 

181  DELVl -XVP(3)  -  XUP<2) 

IF ( ABS( PELU1 )  -  T01 2)153,162,182 

182  N032--N032  +  1 

IF  ( N032  -  KAXNV)154,ie4,200 

183  N032=0 

184  IF (NRES1DI3)  -  NRrlAX )  1 35 •  1 83 , 200 

188  BLSF =3. 

200  IF (BESF 1220,201 ,220 

201  IF(LFLAG(2) 1210,210,202 

202  IF (0RFE(2. 1 )  -  ORFF. ( 1 , 1 )) 203 ,203 , 210 

203  DEL01 =X0P(2)  -  XOP(l) 

IF ( A8S ( DEL01 )  -  T0L2 ) 205 , 204 , 204 

204  NO 21 =  N021  +  1 

IF ( N021  -  MAXNO >208,208,210 

205  N021 =0 

208  IF (NRESID(2)  -  NRMAX)207,207,210 
207  BESF=2. 

GO  TO  220 

210  I F ( LFLAG (1)1212,212,218 

212  IF(LF(_AG<2)  >213,213,215 

213  eES F=3 

GO  TO  220 
215  eESF=2 

GO  TO  220 
210  BESF=1 

220  TIME1=TIME  *  T SC ALE 
IBESF=BESF 

RESID(4)=RESID( IBESF) 

RPSIG(4) =RP5IG( IBESF) 

VRFE(4, 1 ) =VRFE( IBESF , 1 ) 

VRFE ( 4 » 2 ) = VRFE ( I BESF , Z ) 

XPP(4)=XPP< IBESF) 

XVP(4)=XVP< IBESF) 

XAP(4)=XAP( IBESF) 

DO  500  1=1,6 

IPRINT  < I )  =  IPRINT  < I )  +  1 
IFUPRINTU)  .E0.  -1000JG0  TO  254 
IFdPRINT(I)  .LT.  I GAP! I )  )G0  TO  500 

254  IPRINT(I)=0 

GO  TO  (255,285, 270, 275, 27B, 277), I 

255  PRINT  256»TIME1 

256  FQRMAT(/2X,F10.2) 

PRINT  258, (RESID( JR) , JR=1 ,4) ,BESF, <NEX( JR) , JR=1 ,3) , 
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*  (LFLAG<  JR) . JR  =  1 ,3) 

C  WRITE! 8)TIME1 , (RESIDE  JR ) , JR  =  1 , 4 ) .BESF 

258  FORMAT (2X.5E15.7.SI3) 

GO  TO  500 

265  URITE<2>256)TIME1 

UR1 TE(2>  257) ( SORT (RPSIG! JR) ) , JR=1 >4) ,BESF 

257  FORMAT (2Xi 5E15.7) 

C  WRITE (9)TIME1» (SQRT(RPSIG( JR) ) » JR=1 » 4 ) .BESF 

GO  TO  500 

270  WRITE! 3 .256) TIME1 

WRITE (3. 257) ( (ORFE! JR.KK ) > JR=1 .4) ,BESF,KK=1 .2) 

GO  TO  500 

275  WRITE(4.256)TIME1 

WRITE! 4 r 257) (XPP! JR) . JR=1 . 4 ) ,BESF 
GO  TO  500 

27G  WRITE!5,25G)TIME1 

WR I TE  <  5 . 257 ) ( X0 P  <  JR ) , JR= 1 . 4 ) , BESF 

C  WRITE! 1 1 ) TIME1 » (XVP! JR) . JR-1 i4)> BESF 

GO  TO  500 

277  WRITE(6.256)TIME1 

WR1TE(6.257)(XAP(JR),JR=1.4),BESF 

500  CONTINUE 

READ<10)TIME.<W0RD(N),N=1,7Z).<KFG(N>,N=1.12),<TIT(N),N=1,3) 
GO  TO  7 

999  IFINPLQTT  .NE.  'YES')GO  TO  1110 

1110  STOP 
END 
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1  i  me 

F  i  1  ter  1 
Residual 

Filter  2 
Residual 

Filter  3 
Residual 

Best 

F i 1  ter 
Residual 

Best 

F  i  1  ter 
Number 

2  0.90 

1  1  o-i.  8  5  IE  *0  2 

-• 13575?lt+02 

-.  luC8r>32E+G2 

j.l>Cc5p  JcC+02 

.  30OQQO0r+0l 

2  0.90 

5  C  7  3  i  ^  c  t.  +  0  1 

-.  8  194788:  +0  1 

-.12923 39  E +02 

- . 12923 J9 E+ 02 

.  3o(j0000E  +  01 

2  1.00 

£  0  v  j  b  :>  ft  t_  ♦  0  0 

- .  4  30  2493E+  C  1 

-.  6785  02  96  +0 1 

- . 8  785  0  3  9  £  +  0 1 

• JOJdOOOrfOl 

yi .  to 

<^93  )  1  fl  It  ♦  u  1 

-. 3094942L+01 

136  7  4 1 2E +02 

-.  13874  12L  + 02 

.3  JOOOOOE  +  Ol 

21.20 

6  0  2  6  7  9  5  t  +  0  1 

-. 1367421E+02 

-.  2 18  561  BE  +0  2 

-.21  8 561 8  £  +  02 

.  3  0  000  0  0E  +  01 

21.30 

2751791E+01 

-. 3364519E+01 

-. 1217  817E  +  02 

- . 1 217  817E+  02 

.  30Q0000E+01 

21.90 

130502QE  *-01 

-.  4268441E  +01- 

1346946E+02 

-.1305020  E  +01 

.  lOJOOOOE  +  Ol 

21.50 

2  3  8  3  5  9  6E  +01 

-. 1957227E+0 1 

-.  1 1 2  8  9  6  3  E +02 

.  2383b  96E  +  01 

.  300oooo:  +  oi 

21.80 

36320232 +01 

. 1079361 E +01 

77b6589E+0 1 

. ID  79 3  b  IE +  01 

.  2  00  00  0  0E+01 

21.70 

1972979E+Q1 

-. 5 169665 E+ 01 

-.  148  4  365!  +0  2 

-  .  5T69665E  +  01 

.  20  090  DOE +  01 

21.80 

2037006E+Q1 

21.90 - 

1 45827  IE  +01 

5923866E+01 

-.  lbbl 126E  +  02 

- . 5923866E+01 

.  20  000  0  0E  +  01 

-. 579442 5E+0 1 

-.1749429E+02 

-.  5794425E  +  01 

.  200000  OE +  01 

22.00 

113674 9E +00 

- . 431 7745E+0 1 

-.  lo73  432E+02 

-.  43 1 7  7  45E  +  01 

.  20000  0  0E  +  01 

22.10 

22838 54E +00 

-. 4347469E+01 

-.1750  837E+02 

-  .  434  74  t>9E+  01 

.2000000E+01 

22.20 

2  53>  ZZ4.6  E2+ 0  Q _ - .  J*  3-0  J±ZZ3  _ ±  0 1 _ ri.  J.  8.2Z151  EL+ Q_2 -  ,_4  3  fl  4 2 23 £t  ft J}  J _ .2DODOOQE+  0  1 


22.30 

1484098E+Q1  5989348E+01 — 2  1C662  3E +0  2 .  14 840 98E+ 01  .  100000QE+01 


D-3  Range  Measurement 

Adaptive  Exponential  Filter  Output  Example: 
Residuals 

Figure  C.l 


T  ime 

Filter  1 

S  i  gma 

Filter  2 

Siqma 

Fil ter  3 

S  igma 

Best 

F i 1  ter 

Si  gma 

Best 

F i 1  ter 

Number 

20  .  HO 

b4')b727E  +  Ql 

.  9 753  15  OF* Cl 

.  103517  St  +02 

-  .  10  351  7-8E +02 

.  3  0  0000CE  +  01 

20.30 

.  9  70  1 99  5  E  +  0  1 

. 1CE5221E+02 

. 1 C552J1C+02 

.3  flOCOOQt+Gl 

2 1  .  0  3 

61  7  b  92  7E  +  01 

. 9  425253 E  +  C 1 

.  104383  IE +  02 

.  10438  31E  +  02 

.  33000001+01 

21.10 

oQO  73b7:_+01 

.  9  3  3  3  3  0  7  -I  +  0  1 

. 1G71474E+02 

. 1071474L+02 

.  3  0  0  0  0  OOE  +  O  1 

21.20 

o00«g49l+01 

. 97G3G80E+01 

. 1184094E+02 

.  118-.  094E+  02 

.  3C3330GE+01 

21.30 

5fi  <tu  oy  IE  *-01 

. 9406635E+01 

.  US  6485E  +02 

.  1 1  8  o  4  8  5  E  +  0  2 

.3  30330 OE +01 

21.40 

5642553E+01 

. 9141456E+01 

.  1198  417E+02 

• 5642553E+01 

.  1000DQ0E+01 

21 . 50 

547 7907E+G1 

. 833G897E+01 

. 1193687E+02 

. 5477907E+01 

.  13C3G00E+01 

21.50 

535 93906+01 

.  3  5 2  0  9 3 5 £  +  0  1  ... 

.  11L9341E+02 

.  85209'JbE  +01 

.2  30  00  0  01  +  01 

21.70 

5  204297E  +  01 

. 3330404E+01 

. 1194 1 01E+02 

.  83304  J4E  +  01 

.  2  CJ0000E+01 

21 . 80 

50  4763  bE  <-01 

.8185010 E +01 

.  1232567E+02 

.  8185  0  10E+  01 

.  2  0  3  00  0  0E  +  01 

21.90  - 

4683088E+01 

. 8G40837E  +  01 

. 1276582E  +02 

.  8  040  837E  +  01 

.  20  000  0  0E  +  0  1 

22. 00 

470917 faE +01 

. 7337997E+01 

.  13G7379E  +  02 

.  78379  J7E+01 

.230000  3E+01 

22.  10 

4541767E  •■01 

.7645708E+01 

. 1343195E+02 

.  7  645  7  0  8E  +01 

.2  00000  OE +  01 

22.23 

4380454E  *01 

. 7h60678E+01 

. 1382135E+02 

. 746Cb78E+01 

. 2030000E+01 

22.30 

42425bfaE+01 

. 736725 5E+01 

•  1444724E  +0  2 

*  h2425  06E+  01 

. lOOOCOOE+Ol 

D-3  Range  Measurement 

Adaptive  Exponential  Filter  Output  Example: 
Standard  Deviations 

Figure  C.2 
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Position  Time  Best 

'  locity  VRF1  VRF2  ■  VRF3  VRF-Best  Filter 

VRF1  V  R  F  ?  VRF3  VRF-Best  Rumber 

20.-3  0  . 

.  6  <>1  774  9r  *  00  .24993?8E+00  .  1  7  3  0  5  5  61  +  0  0  .17305561+00  .  jJ  DLtL)U00B  +  Dl 
.  2  4  9  3  4  3  B  r.  ♦  0  1  .3062371  1  +  00  .110979  OF +0C  .11097901+00  .30000001  +  01 

20.90 

.  061  77491  +00  .  24993581  +  0  0  .  1697  91  21+00  .16979121+00  .  3  0  00  0  0  01  +  01 

.  2  49  343  0..  :  J 1  .3  0  623711+00  .1U62  13  9F+D0  .  1  0621  356+  00  .  30  0  00  0  01  +  0  1 

p  1  ^  _Q  C 

.  6  6  1  7  7  4  9r ♦ 0  0  .  24993581  +  0  0  .  1  6664761  +  0  0  .  1  6 bu4  76L  +  0 0  .  3  0  000  0  01+  0  1 

.  2  49  3  43".  +  01  .30623711  +  00  .  9  9840441-01  .  9  38  4  0  441-01  .300000  01+01 

21.10 

.  6b 1  7  7 4  91 +  00  .  24993561+  0  0  .  1636  1631  +  00  .  16  361631+  00  .  3  0000001+01 

.  24  934  381  +  01  .30623711+00  .946  2  7091'  -01  .94  82  7u  91-01  .  30000001  +  0  1 

21.20 

.-^61  7  74  91+00  .  2-.993581  +  00  .  1606  9721  +  00  .  16069721+  00  .  3  0  0  G  0  C  0  F  +  0  1 

.24934361+01  .30623711+00  .90143811-01  .90143811-01  . 30000001+01 

21.30 

.55177491+00  .24993581+00  .16787861+00  .157S736E+00  .30000001+01 

.2*9  3-*  -  *  1  +  0  i  .  3  06237  1  1+  00  .85763  891-0  1  .  85763891-  01  .  3  00  00  0  01  +0  1 

21.40  .  _  _  _ 

.66177491+00  .24993581+00  .1551571E+00  .56177491+00  .10000001+01 

.24934381+01  .30623711+00  .B15u3l6E-01  .2+934381+01  .10000001+01 

2  %  r  n 

.  5  61  774  9 E  +  00  .  24993581  +  00  .  1  5252  781  +  00  .  56177491  +00  .  1  0000001  +  01 

.24934381+01  .3062371i+00  .77819721-01  .24934381+01  .10000001+01 

21.60  ... 

.66177491+00  .2+993581+00  .14998616+00  .24993 581+ 00  .20000001+01 

.  249  343  81  +01  .  3  062371  1  +  00  .  7  421  3  711-0  1  .  3  062371  1  +  00  .  2 3 J  0 0 0  Of ♦ 0 1 

21.70 

.56177491+00  .24993581+00  .14752761+00  .2499358E+00  .20000001+01 

.  24934382  +01  .  3  0  b  2  57  1 1  +0  0  ’.  7  0  8  2  7  0  81  -U  1  .3D62371E  +  00  .  200  00  00U-01 

21.80  _ _ _  _  .  _ 

.56177491+00  .24993581+00  .1451487E+00  .2499358E+00  .20000001+01 

.24934381+01  .30623711+00  .67643401-01  .3062371E+00  .20000001+01 

?1 .90 

.5bl7749E+00  .24993581+00  .1428452E+00  .24993581+00  .2000000E+01 

.24934381+01  “  '  .  3  0  6 237  1 1+TJ  0  TbBBVT 6 8E=DT - T3D  6.2  371E  +  DD - V7.  00  0  0  O  01+  0 1 

22.00 _  ...  .  .  _ _ _ _ _ _ _ _ • 

. 561 77491 +00  .2+993581+00  .14061361+00  .24993581+00  .200000GE+01 

.24934381+01  .30623711+00  .6182627E-01  .3062371E+00  .20000001+01 


22.10 

.56177491+00  .24993581+00  .1384506E+00  .2499358E+00  .20000001+01 

-.2  44  3  4  3  81 +TJ1  -25  0  6  2  37  IT  +TTO  V5  9 1T5T6BT=TI  1  .  3  O  62371E+T0  TT00DD00E+D1 


D-3  Range  Measurement 

Adaptive  Exponential  Filter  Output  Example: 
Variance  Reduction  Factors 
Figure  C.3 
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T  i  me 
Filter  1 
Pos i ti on 


2  u  .  BO 

.  2  1  ?  0  C:  j  2  E  *  0  6 
'3  .  90 

?  1  1 9  i  3  7  E  +  0  6 
2  1  .  0  C 

PllnL'u  u  Uu 

21.10 

2  t  1  I  S  7  -rf  +  0  C 

?  1  .  2  0 

FI  lb  49  IE  *36 
21.30 

2  1  1  ■->  b  C  9  t.  +  3  6 
2  1  .  4  0 

2  1 1  4  7  b  6  i  +  0  6 

21.30 

2  11696  /  L  +  06 
2  1  .  b  0 

2113130E+Q6 
2  1.70 

211 229  3E  +  06 

^  2 1  .  9  0 

FI.  90 

211G618E+06 

22.00 

210 180  BE  +06  “ 

22. ID  _  .  .  .  . 
210  8987E+06 

22.20 

210  8180E  +06 

22.30 

210  7J6SE  fOfe 


Fi 1  ter  2 
Pos i Lion 

.2  1  2  0  0  7  1  t  *  3  6 
.  ?  1  1  9  1  7  7  6  +  0  6 
.  2  1  1  8  3  C  L  l  +  0  b 
•2il?-2o;+Go 
.211 692  8  F  +  3 b 
.211  967  21  +  G  6 
.211+819E+06 
.21131816+06 
.  2  1  1  3  1  6 1  ■  +  0  6 
.  2  11 2329  6 +06 
.21114106+06 
.  2  HGfeb  1E+  06 
.2 109542E+06 
V 2 10  9  0  2 8 E  +  c 6 
. 2108220E+06 
.2107410E+G6 


Filter  3 
Posi ti on 

. I 120 C 161+36 
.  c  1 19  2  1  6E  +  0  6 
.1  llBJ6lf_  +  05 
.  2  i  1  7  8  4  ;■  (JC 

.?  116  6  1  0  E  +  ub 
.2119  7992  +  06 
• c 114911C+06 

•  2114074c  +  D6 

•  2113290E+06 
.  2112  42lE  +  06 
.2111  G  9  7 1 +  0  6 
.  211 C-  77 BE +  06 
.22Ct996&E  +  06 
.  21C.9I6  0E+0  6 
.  2108  359E  +  06 
.  210  7561E  +06 


Best 
F i 1  ter 
Posi i i on 

.  2  1200J6E+36 
.  2  1192  :  +  06 

.  .21183912  +  06 
.  2  1 1 7  n  o  4  f .  +  G  6 
.  (.’11661  0L  +  Ub 
.  2119  79  8l  +  06 
.  21  If  7  -j3l  +  Db 
.  2  1 1 3  9  3  7  E  +  3  6 
. 2 1 131 olE +  06 
.  2  1 1  2  3  2  9  L  +  06 
.2111  4  JOE  +  06 
.  2  1 1  u  6  ‘j  1 E  +  0  6 
.  21 098  h2C  +  06 
.  21  090  2  8  E  +  06 
.  2108220E  +  06 
. 21073  l«5E  + 06 


D-3  Range  Measurement 

Adaptive  Exponential  Filter  Outnut  Example: 
Position 
Figure  C.4 
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.  3  jJju00E  +  01 

#  3  0  J  0  0  0  3:.  +  0  1 
.  3  0  JO  JO  Jc +  01 

•  o  0  j  0  G  c  0  f_  +  0  1 

.  30300032+01 
.3  J  J  0  0  0  0:  +  0  1 
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.  2 ooooo or +  oi 

. lOJQOOOE+Ol 
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i  i  me 
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Filter  2 
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8  L 7  0  c  1  5  r +03 

-  .  8  b  1  8  1  7  1  ■:  +  L  3 

2  3.00 

h  7  1  8  3  3  E  +  0  3 

-  .  3  5  7  o  8  0  2  r_  +  0  3 

2  i  ■  j  if 

ri  t  -j  8  7  8  4  f  +03 

-  .  8  >  1  4  H  3  1  -  +  0  3 

2  1.10 

.1  o '  >  r  .  4  0  r  1-03 

-  .  8  u  7  2  4  4  3  E  +  0  3 

21  .  20 

.8  7  3  1  b  3  4 +  LI  3 

-.  3  -t  5  5  9  u  O  E  +  L  3 

'-l  4  -7  0 

(  A  •  J  u 

rt  f  1  5  5  i  a  E  +03 

-.  3  5 4 b  C  5  G  +03 

21.40 

8  5  3 b 2  b 4  E  +  0  3 

-.83399132+03 

21.00 

c  -  4 1  +03 

-.  82741345  +  03 

2  1  .  3  0 

6 2  o 4  2  7  :  L +03 

-  .  H  1  9  4  95  3  E  +  0  3 

21.70 

8  2  u  1 1  3  5  E  +  0  3 

81'+37-»5E  +  03 

21.80 

82220 3 2E +03 

3395PQ3E  +  03 

21.90 

81 9 311  2E +03 

-.  3Q5G276E+03 

22.00 

8135493E+03 

-.79980791+03 

22.10 

60  b  0  73  IE +0  3 

7946744E+03 

22.2  0 

6027116E+03 

-.78959461+03 

22*3  0 

80Q  25  0  0E  +03 

7 8534591 +0  3 

Pest 
F  i  1  ter 
Vel oci ty 

-•SbOGQJ  4  E  +  03 

-  .  853u3  1  8E  +  03 

-  .  34o7UE»8E  +  03 

-  .  340  4  1  3  3E  +  0  3 
-.  83509  V3E  +  03 

-  .  828bl  5bE  <-03 
- . 35282^41+03 
- .  34134u5L  +  03 
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- . 814  3  7  4  5  E  + 03 
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-.7946744E+03 
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-  .  7  b~L  5  1 4  3  E  +  0  3  -.800  25  001+03 
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- .  8  ECO  0  0  4 E  +  0  3 

-.8yJu318E+03 

-  .  8  4 1  7  L.58E+03 

-.84041 8 CE +03 

.83-1  9  4  9 E  +03 
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.B2l27  0  3E  +  O3 

•  8l.5b78?E+G3 

. eo 87 048E+03 
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D-3  Range  Measurement 

Adaptive  Exponential  Filter  Output  Example: 
Vel oci ty 

Figure  C.5 
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.  luooooor+oi 
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.330000  OF +  01 
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. 200000 OE +01 
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6  H  6  b  9  8  4 1  +  "  - 
21.03 
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^1.10 
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21.20 

*>198424:.  +02 

21.30 
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21.40 
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21.-0 
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21.60 
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21.70 

6655427E+02 
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21.90 

602667  3E  +02 

2  2.00 

60062246  +02 

22. 10 

5965141F+02 

22.20 
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22.30 

565171 1E+02 
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.  31  6  •->  1  7  3  r  +  G  2 
.  6  u  -3  3  5  9  8  E  +  0  2 
. 7961 102E+G2 
. 78  2  4705  l  +  02 
. 75 942 4 3E +02 
.  7  3  u  2  3  7  E  +  0  2 
.  7462314E+3-2 
.74293  3  c-:  +  02 
.  7-r*.7E23E  +  02 
. 7360415E+02 
. 7260600E+02 
. 7162965E+02 
. 7090212E+02 
.  7  0  16959E+02 
-.-6  9  4  4  43 3  £t+  0  2 
.  6843514E-+02 
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.  6 1 39 4 65 E +02 

.  8  C  a  4  2  9  6 1  +02 

. oC48876E+02 
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.7917 14 6E +02 
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. 77332G8E+02 
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.  8 1394  6 56  +  02 
.  80  842  9 3E +  02 
.  8  0  4  8 8 7  6  E  +  02 
.  7  995 9  3 2 E +02 
.  79171  46E  +  02 
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D-3  Range  Measurement 

Adaptive  Exponential  Filter  Outout  Example: 
Acceleration 

Figure  C.6 
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STATEMENT  OF  WORK 


During  Task  SER  82-1,  algorithms  required  for  inclusion 
in  the  VAX  11/780  computer  real-time  range  safety  programs  were 
identified.  The  identified  algorithms  were  of  three  types: 
a)  algorithms  being  utilized  in  the  current  (P2457)  software, 
requiring  no  modifications,  b )  algorithms  utilized  in  P  2  4  5  7  , 
requiring  modifications,  and  c)  new  algorithms  not  currently 
included  in  P2457. 

Among  the  new  algorithms  recommended  were  the  following: 

1)  Table  look-up  drag  and  wind  corrected  instantaneous 
impact  predictions  (I  IPs)  with  uncertainty  ellipses. 

2)  Adaptive  exponential  radar  data  filters. 

3)  Utilization  of  telemetered  inertial  navigation 
system  (INS)  and  altimeter  data. 

The  contractor  is  directed  to  develop  the  above  three 
algorithms  and  present  them  in  a  form  conducive  to  fast,  efficient 
incorporation  in  the  VAX  11/780  real-time  range  safety  software. 
The  task  will  proceed  within  the  following  guidelines: 

1)  The  algorithms  will  be  defined  both  in  conventional 
mathematical  and  Fortran  notation. 

2)  Flow  diagrams  of  the  algorithms  will  be  developed 
and  presented. 

3)  The  validity  of  the  algorithms  will  be  demonstrated 
by  running  realistic  test  cases  on  the  AD  CDC 
computer. 


